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ABSTRACT 

We  show  that  the  multigrid  algorithms  of  Brandt  can  be  adapted  to 
solve  linear  complementarity  problems  arising  from  free  boundary  problems. 
The  multigrid  algorithms  are  significantly  faster  than  previous  algorithms. 
Using  the  multigrid  algorithms,  which  are  simple  modifications  of  multigrid 
algorithms  for  equalities,  it  is  possible  to  solve  the  difference  equations 
to  within  truncation  error  using  less  work  than  the  equivalent  of  six 
Gauss-Seidel  sweeps  on  the  finest  grid. 
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SIGNIFICANCE  AND  EXPLANATION 


Several  free  boundary  problems,  (including:  saturated-unsaturated 
flow  through  porous  dams;  elastic-plastic  torsion;  and  cavitating  journal 
bearings)  can  be  formulated  as  linear  complementarity  problems  of  the 
following  type.  Find  a  non-negative  function  u  which  satisfies  pre¬ 
scribed  boundary  conditions  on  a  given  domain  and  which,  furthermore, 
satisfies  a  linear  elliptic  equation  at  each  point  of  the  domain  where  u 
is  greater  than  zero.  We  show  that  the  multigrid  algorithms  of  Brandt, 
(in  which  solutions  are  computed  on  a  series  of  nested  grids)  which  were 
developed  to  solve  boundary  value  problems  for  elliptic  partial  differ¬ 
ential  equations,  can  easily  be  adapted  to  handle  linear  complementarity 
problems.  The  resulting  algorithms  are  significantly  faster  than  previ¬ 
ous  algorithms  in  which  only  one  grid  is  used,  since  the  computation  time 
is  proportional  to  the  number  of  gridpoints  on  the  finest  grid. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
sunmary  lies  with  MRC,  and  not  with  the  authors  of  this  report. 


MULTIGRID  ALGORITHMS  FOR  THE  SOLUTION  OF  LINEAR 
COMPLEMENTARITY  PROBLEMS  ARISING  FROM  FREE  BOUNDARY  PROBLEMS 

»  (n  **  (2) 

Achi  Brandt  '  and  Colin  W.  Cryer  ' 

1.1  INTRODUCTION. 

Several  free  boundary  problems  can  be  reformulated  in  the  form  of  an  (infinite¬ 
dimensional)  LCP  (linear  ccmplementarity  problem):  Given  a  polygonal  domain  0  c  Rn 
with  boundary  30,  and  given  functions  f  and  g,  find  u  (defined  on  0)  such  that 
(in  an  appropriate  weak  sense) 


(a) 

jCu(x) 

1  f  (X), 

X  € 

a  , 

(b) 

u(x) 

>  0, 

x  e 

n  , 

(c) 

u  (x)  t£u(x)  -  f  (x)  1 

-  0, 

X  6 

n  , 

(d) 

u  (x) 

-  g(x). 

X  £ 

312 

where  £  is  a  given  second  order  elliptic  operator.  The  restriction  that  0  is 
polygonal  is  not  essential,  but  suffices  for  our  present  purposes.  We  do  not  write 
(1.1a)  in  the  more  usual  form  -£u(x)  +  f  (x)  ^0  because  we  wish  to  maintain  compati¬ 
bility  with  the  notation  in  previous  papers  by  Brandt. 

Well-known  examples  of  free  boundary  problems  which  can  be  written  in  the  form 
(1.1)  include  porous  flow  through  dams  (a  recent  reference  is  Baiocchi  [1978]),  journal 
bearing  lubrication  (Cryer  [1971a],  Cimatti  [1977])  and  elastic-plastic  torsion  (Cea, 
Glowinski,  and  Nedelec  [1974],  Lanchon  [1974],  Cryer  [1979]).  General  references 

include:  Duvaut  and  Lions  [1976];  Glowinski,  Lions,  and  Tremolieres  [1976],  and 
Cryer  [1977],  Glowinski  [1978];  Cottle,  Giannessi,  and  Lions  [1980];  and  Kinderlehrer 

and  Stampacchia  [1980] . 
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If  f2  is  approximated  by  a  regular  grid  then  the  grid  can  be  divided  into 
N  =  | G |  "interior"  points  G  and  | 3G (  "boundary"  points  8G.  Let  the  grid  size 
be  h.  When  (1.1)  is  approximated  using  finite  differences  on  G,  one  obtains  a 
(finite-dimensional)  LCP: 


(a) 

LU(x) 

1  f  (x), 

x  e 

G 

(b) 

U(x) 

>  o. 

x  e 

G 

(c) 

U(x)  [LU  (x)  -  f  (x)] 

=  0, 

X  € 

G 

(d) 

U(x) 

=  g (x) , 

X  € 

3G 

where  U  (x)  is  an  approximation  to  u(x)  at  the  grid  points  x  e  G  U  3G 
L  is  a  difference  operator  which  approximates  £.  The  coefficients  of  L 
0(h-2). 


(1.2) 


and  where 
are 


2 

By  multiplying  (1.2)  by  h  and  eliminating  the  known  values  of  U(x)  on  3G, 
the  LCP  (1.2)  may  be  written  in  matrix  form 

(a)  AU  <_  b  , 

(b)  1)  >0  ,  (1.3) 

(c)  UT(AU  -  b)  »  0  , 

where  U  is  the  N-vector  of  values  of  U(x)  on  G,  and  A  is  an  N  x  N  matrix 
with  coefficients  which  are  0(1).  Since  we  will  assume  that  A  is  symmetric  and 
negative  definite,  (1.3)  could  be  brought  into  the  canonical  form  for  an  LCP  by  multi¬ 
plying  (1.3a)  by  -1. 

2 

For  example,  if  £  is  the  Laplace  operator  in  R  ,  then  a  possible  choice  for 
L  would  be  the  classical  five-point  difference  operator,  in  which  case  A  would  be 
a  matrix  with  diagonal  elements  -4  and  off-diagonal  elements  either  0  or  1. 

The  general  structure  of  a  finite-dimensional  LCP  is  that  we  have  a  pair  of  vector 
inequalities  together  with  the  complementarity  condi tion  which  states  that  at  every 
point  at  least  one  of  the  inequalities  must  in  fact  be  an  equality. 

There  is  an  extensive  literature  on  the  (finite-dimensional)  LCP  (see  Balinski 


and  Cottle  [1978]).  In  particular,  if  A  is  negative  definite,  as  we  assume,  then 
there  exists  a  unique  solution  to  (1.2)  and  (1.3). 


Since  the  LCP  (1.3)  arises  from  a  free  boundary  problem,  the  matrix  A  has 
special  properties  which  make  it  possible  to  use  specialized  algorithms  which  are 
particularly  efficient.  Such  algorithms  include  projected  SOR  (Cryer  [1971],  Glowinski 
[1971])  the  method  of  Cottle  and  Sacher  [1977],  and  the  modified  block  SOR  (MBSOR) 
method  of  Cottle,  Golub,  and  Sacher  [1978] j  Cryer  [1979a]  summarizes  these  algorithms 
and  Cottle  [1974]  gives  numerical  comparisons  between  them. 

Recently,  it  has  been  found  (Brandt  [1977],  Brandt  and  Dinar  [1979])  that  multigrid 
algorithms  are  an  effective  tool  for  solving  linear  equations  of  the  form 

AX  -  b  .  (1.4) 

The  basic  idea  of  these  multigrid  algorithms  is  to  compute  on  a  sequence  of  nested 
grids.  The  computation  proceeds  on  a  particular  grid  until  the  error  becomes  smooth 
and  the  rate  of  convergence  slows,  at  which  point  the  computation  is  transferred  to  a 
coarser  grid.  When  the  error  has  been  reduced  on  tht  coarser  grid,  the  solution  on  the 
finer  grid  is  corrected  using  interpolated  values  from  the  coarser  grid. 

In  this  paper,  we  show  how  the  multigrid  algorithms  FAS  and  FMG  of  Brandt  can  be 
modified  to  solve  the  LCP  (1.3).  We  find  that  the  modified  multigrid  algorithms  are 
substantially  faster  than  previous  algorithms.  Indeed,  with  only  minor  modifications, 
the  standard  multigrid  programs  solve  the  LCP  with  essentially  the  same  efficiency  as 
is  attained  for  linear  equations. 

The  paper  is  organized  as  follows.  In  Section  2,  we  describe  PFAS,  the  projected 
full  approximation  scheme  for  solving  (1.3);  PFAS  combines  the  concepts  of  multigrid 
algorithms  with  those  of  projected  SOR.  In  Section  3,  we  discuss  the  implementation 
of  PFAS,  and  in  Section  4,  we  give  numerical  results  obtained  using  PFAS.  In  Section  5, 
we  discuss  alternative  implementations  of  PFAS,  the  last  of  which  leads  to  substantially 
improved  convergence  (we  also  include  several  unsuccessful  implementations  because  they 
are  instructive). 

In  Section  6,  we  describe  results  obtained  using  PFMG,  the  projected  full  multi¬ 
grid  algorithm  for  solving  (1.3).  The  basic  idea  of  PFMG  is  to  compute  the  initial 


approximation  on  each  grid  by  interpolating  an  accurate  solution  on  the  next  coarser 
grid.  Using  PFMG  we  are  able  to  solve  the  1<CP  to  within  truncation  error  using  less 
work  than  the  equivalent  of  six  Gauss-Seidel  sweeps  on  the  finest  grid. 

Our  results  are  summarized  in  Section  7  and  some  possible  extensions  are 
mentioned.  Finally/  listings  of  the  programs  are  given  in  the  appendices. 


THE  PROJECTED  GAUSS-SEIDEL  ALGORITHM 


It  is  possible  to  solve  the  LCP’ s  (2.2)  and  (2.3)  using  the  projected  Gauss-Seidel 

algori thm  which  we  now  describe, 
k  0 

Let  u  '  (x)  be  an  approximate  solution  of  (2.2)  and  (2.3).  We  compute  recur- 

k  1  k  2  k  s-1 

sively  a  sequence  of  approximations  u  '  (x) ,  u  '  (x) , . . . ,  as  follows.  I^t  u  '  (x) 

k  s 

be  given.  From  (2. 2d),  the  boundary  values  of  u  '  (x)  are  equal  to  g(x).  The 
k  s 

interior  values  of  u  '  (x) ,  which  together  comprise  the  vector 


uk'S  =  (uk'S  :  1  <  j  <  =  {uk'S(xk>  ••  1  <_  j  <_  N^}  ,  (2.5) 

are  obtained,  point  by  point,  by  first  applying  the  classical  Gauss-Seidel  method  to 
(2.3)  to  obtain 


uk,s-i 

3 


k,s-l  ,.k  r  k  k,s 

3  3  z<j  3*  * 

k, s-1  -k, s  .  k 

u .  +  r .  /a . .  ,  say , 

3  3  33 


I 

*2 


k  k,s-l 
ajiuJt 


(2.6) 


and  then  projecting: 

uk,S  «  max{0,uk,S  .  (2.7) 

k  s  k  s-1 

The  process  of  applying  (2.6)  and  (2.7)  for  1  <  j  <  ^  to  obtain  u  '  from  u 

will  be  called  a  Gk  projected  Gauss-Seidel  sweep,  or  a  projected  sweep.  The 

~k  s 

quantities  r_.  '  will  be  called  the  dynamic  residuals. 

It  is  known  (Cryer  [1971],  Glowinski  [1971])  that  uk’s  -*■  U*  as  s  ■+  ~. 

When  implementing  the  projected  Gauss-Seidel  method  only  the  latest  values  of  the 
solution  are  stored.  We  will,  therefore,  often  suppress  the  iteration  counter  s  and 
denote  one  projected  Gauss-Seidel  sweep  applied  to  (2.2)  and  (2.3)  by 


k 

u 


*■  Projected  Gauss-Seidel 


.  k  Tk  k 
[u  :  L  ,F  ) 


(2.8) 


Similarly , 


_  k  k,s  k,s-l 
Vu  =  u  -  u 


(2.9) 


will  denote  the  difference  between  the  latest  approximation  u  and  its  predecessor, 


denotes  the  previous  difference. 

ERROR  ESTIMATES  FOR  THE  PROJECTED  GAUSS-SEIDEL  ALGORITHM 

When  implementing  the  projected  Gauss-Seidel  algorithm  as  part  of  a  multigrid 
process,  it  is  important  to  be  able  to  estimate  the  error.  In  order  to  do  so,  we  note 
that  since,  by  assumption,  -A  is  symmetric  and  positive  definite,  there  exists  a 
coercitivity  constant  >  0  such  that 

wT(-Ak)w  >  o^w  ,  (2.11) 

\ 

for  all  w  e  r  . 

Lemma  2.1 

k  k 

Let  u  be  the  solution  of  the  LCP  (2.3),  and  let  u  >  0  be  am  approximate 


solution.  Let 


rk  =  (rk)  -  bk  -  Akuk  , 

(2.12) 

k 

and  r+ 

,  k  \ 

-  'V' 

where 

rk  if  uk  >  0  , 

k  3  3 

r+j  ’  k  k 

min{0,r  . },  if  u  .  •  0  . 

J  J 

(2.13) 

Then 

(uk 

-  uk)T(-Ak)(Uk  -  uk)  <  (Uk  -  uk)T(-rk)  . 

(2.14) 

Hence , 

II uk  -  ukii2<  a;1iirkii2 . 

(2.15) 

Proof : 

With 

+ 

defined 

as  above,  we  see  that  u  satisfies  the  LCP: 

(a) 

k  k  .  k  k 

A  u  £  b  -  r  , 

(b) 

uk  >_  0  , 

(2.16) 

(c) 

(uk)T(Akuk  -  bk  ♦  rk)  =  0  . 

k  T 

Following  Falk  [1974]  we  multiply  (2.3a)  by  the  non-negative  vector  (u  )  and  use 
the  complementarity  condition  (2.3c)  to  obtain 


I 


k  „k.T  kk  .  k  k.T  k 
(u  -  U  )  A  U  <(u  -  U  )  b 


(*) 


Similarly,  multiplying  (2.16a)  by  (UJc)T  we  obtain 

„Jc  k.Tk  k  „Jc  k.T,.k  k^ 

(u  -  u  )  A  u  _<  (tr  -  u  )  (b  -  r  )  . 

Adding  (*)  and  (**)  and  combining  terms  we  obtain  (2.14)  and  hence  (2.15). 

Lenina  2.2. 


(**> 

□ 


k  k 

Let  U  be  the  solution  of  the  LCP  (2.3),  and  let  u  ^  0  be  an  approximate  solu- 

)t 

tion  obtained  after  one  or  more  G  projected  sweeps.  Let 


(Dk  -  Lk  -  Pk) 


(2.17) 


k  k  k 

where  D  is  diagonal,  and  L  and  P  are  strictly  lower  and  upper  triangular 
matrices,  respectively. 


Then  u  satisfies  the  LCP 


k  k  k  _k  k 
A  u  <  b  -  P  ?u 


k 

u  >  0  , 


(2.18) 


,  k.T,  k  k  .k  _k„  k. 

(u  )  (A  u  -  b  +  P  7u  )  =  0 


Hence, 


lluk  -  uk||2<  a-1||Pk||2||Vuk]|2 


(2.19) 


Proof;  Consider  the  projected  Gauss-Seidel  method  defined  by  (2.6)  and  (2.7).  For 

k  -»k  •  s  k  s 

each  point  x.  we  first  compute  the  dynamic  residual  r.  .  The  new  value  of  u.' 

k 

is  chosen  so  as  to  reduce  the  residual.  Denote  the  residual  at  the  point  x 

<*k  s 

iirmediately  after  step  (2.7)  by  r/  ,  80  that 


*k,s  -k,s  k  ,  k,s  k,s-l. 
r  =  r .  -  a  .  .  (u .  -  u.  )  . 

j  3  33  3  3 


(2.20) 


two  possibilities: 


k ,  s 

and 

-k ,  s 

33 

u, 

3 

>  0 

r.' 

3 

=  0 

k,s 

V 

=  0 

and 

*k,s 

j 

1  0 

'■'k  » k 

Thus,  dropping  the  superscript  s,  and  setting  r  ={r.  :  1  £  j  <  N.  }, 

3  k 


-8- 


uk  >  0  , 


rk  >  0  , 


(2.21) 


/  fc.T-k 

(u  )  r  =  0  . 


Let 


k  k  k  k 
r  =  b  -  A  u 


It  is  readily  seen  from  (2.17)  that 


k  .k  _k ,  k,s  k,s-l. 
r  =  r  +  P  (u  -u  ), 


(2.22) 


*k  _k„  k 
=  r  +  P  Vu  . 


Combining  (2.21)  and  (2.22)  we  obtain  (2.18).  Comparing  (2.16)  and  (2.18)  we  see 

k 


that  the  arguments  which  led  to  (2.15)  from  (2.16)  may  be  applied  to  (2.18),  with  r 


k  k 

replaced  by  -P  Vu  ,  to  obtain  (2.19). 


□ 


As  Leninas  2.1  and  2.2  show,  we  can  estimate  the  error  in  an  approximate  solution 


k  k  k  ^ 

u  in  terms  of  the  residual  r  or  the  difference  Vu  ;  we  will  usually  use  Vu  to 


estimate  the  error,  since  this  quantity  is  readily  available  during  a  G  projected 
sweep . 


k 


Remark.  The  reader  may  wonder  why  we  bothered  to  introduce  r+  in  Lemma  2.1,  since 


k  k 

(2.15)  holds  with  r+  replaced  by  r  .  The  reason  is  that  for  the  LCP  (2.3)  there 


k  k  k 

may  be  large  positive  residuals  at  points  x_.  where  U  (x.)  =  0,  but  this  does  not 


mean  that  the  error  is  large.  O 

In  multigrid  algorithms  it  is  necessary  to  compare  norms  on  different  grids.  We, 
therefore,  wish  to  introduce  a  norm  which  is  not  grid  dependent.  To  do  so,  we  proceed 
as  follows. 


We  first  note  that,  to  a  good  approximation ,  the  coercivity  constant  a^_  for  -A 


satisfies 


ak  =  ah  , 


where  a  is  the  smallest  eigenvalue  of  £. 


Next,  assume  that  the  approximate  grid  function  u  has  been  extended  to  a 


function  u  (x)  on  Q  approximating  the  solution  u(x)  of  (1.1).  Then 


-9- 


r 


i 


|u(x)  -  u  (x) 


2.0 


|u(x)  -  u~(x)  dx 


r  .  n |  k  k 

l.  "klWj  -  Uj 


k|2|4 


1  j-1 


The  norms 
formula. 


“  - ht 

a  K 

are  essentially  independent  of  k; 
2.  Thus  a  measure  for  the  error  || 


for  example,  for  the  five-point 

u(x)  -uk(x)  |j  _  is  provided  by 
2tu 

,  ,  (2.23) 


and  this  norm  will  be  used  in  the  computations. 

PFAS  (PROJECTED  FULL  APPROXIMATION  SCHEME) . 

PFAS  (Projected  Full  Approximation  Scheme)  obtains  an  approximation  u  to  the 

solution  t/*  on  the  finest  grid  GM  by  recursively  generating  a  sequence  of  approxi- 

-k  k 

mations  u  on  the  grids  G  . 

Each  u  is  an  approximate  solution  to  an  LCP  of  the  form  (2.2)  with  F  replaced 

-.)*  -]( 
by  a  function  F  which  is  defined  later.  In  general,  F  is  different  from  F  so 

b  **M  M  «.M 

that  u  is  not  an  approximation  to  u  .  However,  F  =  F  and  so  u  is  an 
M 

approximation  to  U  . 

We  begin  by  initializing  uM  to  some  suitable  value.  For  example,  we  might  set 


uM  (x)  =  g  (x) ,  on  3GM  , 
uM (x)  =0  in  GM  . 


(2.24) 


We  also  set 

[|7uM||g=  lo30,  eM  =  €  ,  (2.25) 

(where  e  is  the  desired  accuracy  on  the  finest  grid,  and  where  the  astronomical  number 
103^  ensures  that  at  least  two  GM  projected  sweeps  are  carried  out). 


-10- 


-M  ,  .  M  .  .  .  M 

F  (x)  -  F  (x),  for  X  «  G  , 

and  (2.26) 

i?1  (x)  »  (Ax),  for  x  e  GM  . 

M 

He  now  make  a  number  of  G  projected  sweeps, 

uM  Projected  Gauss-Seidel  [uM:LM,FM]  .  (2.27) 

After  each  sweep  we  test  whether 


'  M 

If  so,  the  accuracy  criterion  is  satisfied,  and  we  accept  u  as  an  accurate  approxima¬ 
tion  to  UM  *  5”  on  GM. 

It  is  known  that  Gauss-Seidel  iteration  is  a  smoothing  process:  the  error 
~M  -  m 

U  (x)  -  u  (x)  becomes  smoother  as  the  number  of  sweeps  increases,  while,  at  the  same 
time,  the  rate  of  convergence  slows  down.  We,  therefore,  carry  out  only  a  few  gM 
projected  sweeps,  stopping  when  either  (2.28)  is  satisfied  or 

||vGmIIg1  ni|7S“la||G  •  (2.29) 

Here,  n  is  a  fixed  parameter;  in  our  work  we  have  taken  n  “  -5. 

Suppose  that  (2.28)  is  not  satisfied  but  that  (2.29)  i3  satisfied.  This  means  on 

the  one  hand  that  the  accuracy  of  u  must  be  improved  and  on  the  other  hand  that  it 

M  M 

is  inefficient  to  continue  iterating  on  G  .  The  slow  rate  of  convergence  on  G 
indicates  that  the  error  is  smooth,  so  that  the  error  can  be  represented  satisfactorily 

M— 1  M-l 

to  the  next  coarsest  grid,  G  .We  therefore  move  to  G 

-71  M  -M 

Since  IT  (x)  satisfies  (2.2),  with  k  =  M  and  F  =  F  ,  the  error 

Ax)  -  Ax)  -  Ax)  ,  (2.30) 

satisfies  the  LCP 

lV(x)  <  Ax),  on  GM  , 

Ax)  +  uM(x)  >  0,  on  GM  , 

(2.31) 

(Ax)  +  Ax)(  (LMAx)  -  Ax>)  =  0,  on  GM  , 

V1^  (x)  =  0  ,  on  3GM  , 
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where  the  residual  rM  is  given  by 

— H ,  .  -M  ,  .  M-M  .  M  _ . 

r  (x)  «  F  <x)  -  L  u  (x),  x  e  G  .  (2.32) 

M 

As  already  observed,  V  (x)  is  a  smooth  function  and  may,  therefore,  be  accurately 
represented  on  G1*-1.  Furthermore,  con$>aring  (2.31)  and  (1.1)  we  see  that  V^x)  is 
an  approximation  to  the  continuous  solution  v (x)  of  the  LCP 


Iv  (x)  <  rM  (x) ,  x  e  S  , 


v  (x)  +  uM(x)  >  0, 


[v(x)  +  uM(x)l  (JCv(x)  -  rM(x))  =  0 


x  e  fi  , 


(2.33) 


x  e  0  , 
v(x)  =  0,  on  3!)  , 

(where,  by  abuse  of  notation,  rM(x)  and  uM  (x)  are  defined  on  (2  by  appropriate 
interpolation  between  the  values  of  r*1  and  uM  on  the  gridpoints  of  GM) .  Thus ,  a 
good  approximation  to  Ax)  may  be  obtained  by  solving  the  finite  difference  approxi¬ 
mation  to  (2.33)  on  GM  1.  That  is,  vAx)  is  closely  approximated  on  GM  1  by  the 
solution  A  1 (x)  of  the  LCP, 

L  W  (x)  <  S„  r  (x),  on  G  , 


(a) 

(b) 

(c) 

(d) 


AI(x)  +  l”“1uM(x)  >0,  on  GM_1  , 


(2.34) 


(A1(x)  +  A1Sm(x)hlm_1wm_1(x) 


.M-l-M ,  „  M-l 

S  r  (x)]  =  0,  on  G 

M 


A  1  (x)  =  0,  on  3G1 


M-l 


Here  ^  and  s‘‘  *  are  operators  taking  grid  functions  on  G"  into  grid  functions 


on  G 


M-l 


(As  an  aid  in  memorization,  note  that  in  1^  the  subscript  M  and 

M 


superscript  M  "cancel”.) 

M-l  M-l 

The  operators  I  and  S  can  be  defined  in  many  ways.  One  choice  is  to 
M  M 

.M-l 


choose  both  IM  and  ^  to  be  the  injection  operator: 
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(2.35) 


t  .M-l  ,  ,  .  .  ,  M-l 

w(x)  =  w(x),  x  c  G 

M 

.  M-l  M-l 

Other  choices  for  I  and  S  will  be  discussed  later. 

M  M 

If  we  were  solving  a  linear  boundary  value  problem  then  condition  (2.34b)  would  not 

M-l  M-l 

apply  and  it  would  be  most  efficient  to  solve  for  the  correction  w  on  G  Since 

we  are  solving  inequalities  the  problem  is  nonlinear  and  it  is  necessary  to  solve  for  a 
'full  approximation’  uM  *  on  GM 


Setting 


UM-1(x)  =  /^(X)  + 

M-l-M,  . 

!m  u  (X)  , 

(2.36) 

it  follows  that 

-M-l 

U  (x)  satisfies  the  LCP 

(a) 

L^’hx) 

1  FM-1(X), 

M-l 

in  G  , 

(b) 

SM_1 (x) 

>  o, 

M-l 

in  G  , 

(2.37) 

(c) 

UM-1(x)  (t,M'1UM“1(x)  -  FM_1(X)1 

-  0, 

in  G^1  , 

(d) 

uM-1(x) 

»  g  (x) ,  on 

3GM“1  , 

where 

=M-1,  ^ 

F  (x)  » 

^(x)  +lm'1iJJ’1Gm(x)  =  s”'1 

[F  (x)  -  L  u 

,  ,,  ,  tM-1  M-l-M  ,  , 
(x)]  +L  I  u  (x)  . 

M 

(2.38) 

Finally,  we  set 

e"'1  «  6  I|VuM||g  ,  (2.39) 

and 

-M-l  M-l-M  ,,  .  . 

u  =  XM  u  ,  (2.40) 

where  6  is  a  constant;  in  our  computations  6  has  been  set  equal  to  .15. 

-M  M  -M 

To  recapitulate,  starting  with  initial  values  of  u  ,  e  ,  and  F  ,  we  first 

carry  out  GM  projected  sweeps  until  convergence  slows  down.  We  then  introduce  a 

subsidiary  problem  on  G11  1  with  known  FM  1  and  eM  1  and  initial  approximation 
— M— 1 

u  .  The  process  can  be  repeated,  so  that  at  any  one  stage  of  the  computation  we  have 

a  sequence  of  grid  approximations  u^u1,1  _ .u*4  (approximating  UM,0^  *,...,U* 

M  M-l  k— 1  “M  — M— 1  — k— 1 

respectively),  tolerances  e  ,e  ,  ...,e  ,  and  right  hand  sides  F  ,F  ,...,F 
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* 


* 


■ilKSb 


In  the  general 

case,  U*  is  the  solution  of  the 

LCP 

(a) 

LkUk(x)  <Fk(x), 

In 

Gk  , 

(b) 

uNx)  >  0, 

in 

Gk  , 

(c) 

uNx)  (LkUk(x)  -  Fk(x))  -  0, 

in 

Gk, 

(d) 

uNx)  =  g(x) 

on 

,-k 

3G  ; 

or  equivalently. 

(a) 

Ak5k  <  bk  . 

(b) 

uk  >  0  , 

(c) 

(Uk)T(AkUk  -  bk)  =  0 

(2.41) 


(2.42) 


This  LCP  is  solved  approximately  using  G  projected  sweeps  until  the  latest  approxi- 
mation  u  satisfies  either 


(2.43) 


l’uk|L>n||v^ 


(2.44) 


k-1 


'G-""  'old  G  * 

If  (2.44)  holds  but  (2.43)  does  not  then  a  new  problem  on  G"  -  is  defined  by 
setting 

-k-1  Jc-1  -k  Tk-k  k-lk-l-k 

F  =  (F  -  L  u  ]  +  L  u  ,  (2.45) 


e 


k-1 


-k-1 

u 


(2.46) 

(2.47) 


-k-1  Jc-1  Tk-l-k 

U  =  W  +  Ik  u  , 

,,k  nk  -k 

V  =  U  -  U  » 

k-1  k  k-1 

where  W  is  an  approximation  to  V  on  G  .  Unless  otherwise  indicated, 

k-1  k-1 

and  S.  will  be  taken  to  be  the  injection  operator  Inj. 
k  k 

-k-1 

At  some  stage  the  latest  approximation  u  must  satisfy,  (2.43) 


(2.48) 

(2.49) 


Iuk'1| 


k-1 


G  - 


(2.50) 
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if  for  no  other  reason  than  that  when  k-1  »  1  we  cannot  introduce  any  more  sub 


sidiary  problems  and  must  iterate  until  (2.50)  is  satisfied.  Having  found  an  approxi- 

-k-1  x 

station  u  of  sufficient  accuracy,  we  return  to  G  .  To  do  so,  we  first  determine 

X-i  X-l 

an  approximation  w  to  w  from  (2.48)  namely 

k-1  -k-1  Jc-l-k 

w  -  u  -  lk  u  .  (2.51) 

k  k-1 

Next,  let  IJc  ^  be  an  interpolation  operator  taking  grid  functions  on  G  into 

k  k 

grid  functions  on  G  .  A  possible  choice  for  I  ^  is  the  linear  interpolation  operator 

^  defined  as  follows.  If  P^,  P^,  P^,  and  P4  are  the  corners  of  a  square  in 
k-1 

G  (see  Figure  2.1)  then 


,k  k-1  . 

Lk-lW  (V 


wk_1(Pi),  1  <  i  <  4  , 

(w1C_1{P1)  +  wk_1(P2))/2,  i  -  5 
(wk"1(P1)  +  wk_1(P4))/2.  i  -  6 


(  1  wk_1(P  ))/4,  i  -  7  . 
i-1 


(Other  choices  for  ^  will  be  discussed  later.) 


is  an  approximation  to 


is  an  approximation  to  v  ,  and,  noting  <2.49), 

-k  -k  „k  k-1  ... 

u  «  u  +  I^_^w  ,  (2.54) 

is  an  in?) roved  approximation  to  U  .  However,  because  of  the  nonnegativity  constraint 

upon  U^,  we  allow  somewhat  greater  generality  and  replace  u*  as  follows: 

-k  ,~k  -k,  ,-k  ,  k  k-1  -k.  ... 

u  -*-V>(u  ;u  )  “  V>(u  +  ;u  )  •  (2.55) 

Initially  we  set 

,-k  -k.  ~k  ... 

^>(u  ;u  )  »  u  ,  (2.56) 

but  other  choices  will  be  considered  later. 

PFAS  is  described  by  (2.24)  through  (2.56).  A  flowchart  is  given  in  Figure  3.1, 

and  the  implementation  is  discussed  in  Section  3.  If  the  algorithm  converges,  we  will 

_M 

eventually  obtain  an  approximation  u  satisfying  the  required  accuracy  condition 
(2.28)  and  the  algorithm  will  terminate. 
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3.  IMPLEMENTATION  OF  PFAS. 

The  flowchart  for  PFAS  is  given  in  Figure  3.1.  PFAS  has  been  implemented  as  a 

2 

FORTRAN  subroutine  for  the  case  when  £2  is  a  rectangle  in  R  ,  £  is  the  Laplacian 

k-1  k-1  k 

operator,  1^  and  are  injections  (equation  (2.35)),  and  1^  ^  is  linear 

interpolation  (equation  (2.52)).  The  subroutine  PFAS,  which  is  listed  in  Appendix  A 

as  part  of  the  program  for  solving  the  porous  flow  free  boundary  problem  described  in 

Section  4,  is  a  straightforward  modification  of  an  earlier  program,  FAS  Cycle  C,  of 

Brandt.  In  the  subroutine  PFAS  most  of  the  computations  are  performed  by  auxiliary 

subroutines,  and  the  flowchart  shows  the  role  played  by  these  auxiliary  subroutines. 

One  reason  for  giving  a  listing  of  PFAS  is  so  that  the  reader  can  appreciate  how 

easy  it  is  to  implement  PFAS.  It  may  also  be  remarked  that  many  other  interesting  free 

boundary  problems  (for  example,  elastic-plastic  torsion  problems  and  cavitating  journal 

bearing  problems)  are  formulated  on  simple  polygonal  regions,  and  the  program  given 

here  could  easily  be  modified  to  handle  these  problems. 


The  following  coiments  arise : 

-k 

1.  In  PFAS,  the  LCP  for  D  is  solved  in  the  form  (2.42)  rather  than  (2.41), 

-k  k  -k  2-k 

but  the  values  of  u  on  3G  are  also  stored.  Thus,  b  =  h, F  is  stored  instead 

k 

-k  k  k-1 

of  F  .  In  going  from  G  to  G  we  have,  from  (2.45),  since  ^  =  2hk» 

-k-1  2  -k-1 

b  =  hk.lF  , 

=  h^_i  (S^'*  [Fk  -  Lkuk]  +  L-lTuN  . 


k-1  -2rr-k  k-k. 


k-1  k-1 -k. 


=  hk  [b  -  A  u  ]  +  L  Ik  u  )  ,  (3.1 

.  k-1  -k  k-k  k-1  k-l-k 

=  4Sk  (b  -  A  u  ]  +  A  1^  u  . 

k  k 

2.  A  G  -work-unit  is  the  work  required  for  one  G  projected  sweep.  The  work 

k  “ n ( M-k )  M 

for  one  G  projected  sweep  is  approximately  2  'g  -work-units,  and  PFAS  keeps 

M 

track  of  the  total  number  of  G  -work-units,  WU.  When  no  confusion  is  possible  we 


write  "work  unit"  instead  of  "G  -work-unit". 
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ft,... 


factor  y,  which  is  defined  by 


2=  lira  f||PuM||  }1/WU  (3.2) 

WU-»~ 

4.  All  the  numerical  computations  were  performed  on  the  Univac  1180  at  the 
University  of  Wisconsin-Madison.  The  programs  were  written  in  ASCII  FORTRAN  and  com¬ 
piled  and  executed  using  full  optimization. 

The  Univac  1180  single-precision  arithmetic  has  approximately  eight  decimals. 

The  residuals  usually  decrease  quite  rapidly  at  the  beginning  of  a  computation  so  the 

round-off  threshhold  is  quickly  reached.  For  example,  for  the  problem  considered  in 

Section  4  with  M  =  5,  ||uM||  is  about  2x10^  and  the  single  precision  algorithm  went 

G 

into  a  loop  when  ||VuM||  reached  5x10  6  after  a  mere  50  work  units. 

G 

In  the  numerical  experiments  we  were  particularly  interested  in  measuring  the 

asymptotic  convergence  factor  y.  TO  eliminate  round-off  effects,  all  the  computations 

reported  on  here  used  double  precision  arithmetic.  Of  course,  this  is  not  normally 

necessary.  Furthermore,  even  if  very  accurate  solutions  of  the  discrete  problem  (2.2) 

-M 

were  required,  it  would  suffice  to  store  u'  in  double  precision  and  all  other 
quantities  in  single  precision. 

The  execution  times  quoted  are  those  provided  by  the  Univac  1180  Exec.  System. 

As  is  often  the  case  on  timesharing  systems,  the  times  are  only  reproducible  to  within 
about  10%. 

Because  of  its  word  length,  the  UNIVAC  1180  can  only  directly  access  64k  words 
of  storage.  When  M  ^  7  more  than  64k  words  of  storage  are  needed  by  PFAS  and  there 
is  a  significant  degradation  in  performance. 

5.  To  measure  y  the  iterations  were  continued  for  the  first  100  work  units, 
unless  the  residuals  vanished  before.  In  practice  one  usually  iterates  only  for  about 
30  work  units. 

We  also  used  several  values  of  M  in  order  to  measure  the  dependence  of  y 
upon  M. 
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I 


Part  of  the  output  of  a  typical  computation  using  PFAS  is  shown  in  Figure  3.2. 

After  each  G  projected  sweep,  the  values  of  the  level  k,  the  residual  norm  ||Vu  ||  , 

G 

and  the  number  of  work  units  WU  are  printed  out. 

The  computations  starting  at  a  level  M/level  (M-l)  junction  and  continuing 
until  the  next  level  fr^/level  (M-  1)  junction  are  called  a  cycle  (see  Figure  3.2). 

For  the  cycle  shown  in  Figure  3.2,  ||VuM||  decreased  from  .293  10~9  to  .110  lo’9 

G 

with  the  expenditure  of  (99.039-94.400)  =  4.639  work  units. 

While  minor  variations  do  arise,  a  cycle  often  consists  of  a  sequence  of  2  sweeps 
at  each  of  levels  M  -  1,M  -  2,...,1,  followed  by  2  sweeps  at  each  of  levels 
2,...,M  -  1,  terminating  with  2  or  3  sweeps  at  level  M.  If  this  pattern  is  followed 
with  3  sweeps  at  level  M  then  the  average  number  of  work  units  per  cycle  is 

3  +  4  [2_n  +  2-2n  +  ...]  =  3  +  4/(2n  -  1)  ,  (3.3) 

and  the  average  number  of  work  units  per  GM  projected  sweep  is  1  +  4/ (3(2°  -  1)). 

Of  course,  very  irregular  patterns  are  observed  when  the  round-off  threshhold  is 
reached. 

6.  As  can  be  seen  from  Figure  3.2,  f( VuM f (  decreases  steadily  but  not  very  regu- 
larly,  in  part  because  of  slight  variations  in  the  number  of  sweeps  at  each  level.  To 
evaluate  the  algorithm,  we  have  used  two  quantities: 


rf  =  H^finalUc 


=  [life 


the  value  of  || 
cycle  before  100  work  units, 
-M 


-M  n 

Vu  I)  at  the  end  of  the  last  complete  (3.4) 
G 


final  Mg'  II’ ’“initial  K 


1/[WU  .-WU.  ...  ,] 

final  initial 


(3.5) 


where  II  ^uj*nitial  ^G  va^-ue  °£  ll7uM|IG  after  the  first  GM  sweep,  is  an 

O 

estimate  for  the  asymptotic  convergence  factor  y. 

For  example,  for  the  data  in  Figure  3.2,  the  value  of  II  Vuin|tial  llG  (which  is 

not  shown  in  Figure  3.2)  was  4.95  and,  of  course,  WU.  .  .  ,  =  1.  Thus, 

initial 


LEVEL  5  RESIDUAL  NORM=  .755-010  WORK=  91.400 

LEVEL  6  RESIDUAL  NORM=  .126-008  WORK=  92.400 

LEVEL  6  RESIDUAL  NORM-  .515-009  WORK=  93.400 

LEVEL  6  RESIDUAL  NORM=  .293-009  WORX=  94.400 

********************END  OF  CYCLE** ************** ****MU  =  .7771 

LEVEL  5  RESIDUAL  NORM=  .196-009  WORK=  94.650 

LEVEL  5  RESIDUAL  NORM=  .133-009  WORK=  94.900 

LEVEL  4  RESIDUAL  NORM=  .879-010  WORK-  94.963 

LEVEL  4  RESIDUAL  NORM=  .613-010  WORK=  95.025 

LEVEL  3  RESIDUAL  NORM-  .385-010  WORK=  95.041 

LEVEL  3  RESIDUAL  NORM=  .257-010  WORK=  95.057 

LEVEL  2  RESIDUAL  NORM-  .133-010  WORK=  95.061 

LEVEL  2  RESIDUAL  NORM-  .717-011  WORK=  95.064 

LEVEL  1  RESIDUAL  NORM=  .243-011  WORK55  95.065 

LEVEL  1  RESIDUAL  NORM-  .447-012  WORK=  95.066 

LEVEL  2  RESIDUAL  NORM555  .303-011  WORK55  95.070 

LEVEL  3  RESIDUAL  NORM55  .189-010  WORK=  95.086 

LEVEL  3  RESIDUAL  NORM=  .714-011  WORK=  95.102 

LEVEL  4  RESIDUAL  NORM=  .686-010  WORK=  95.164 

LEVEL  4  RESIDUAL  NORM55  .255-010  WORK=  95.227 

LEVEL  4  RESIDUAL  NORM=  .138-010  WORK=  95.289 

LEVEL  5  RESIDUAL  NORM=  .151-009  WORK=  95.539 

LEVEL  5  RESIDUAL  NORM-  .5 34-010  UORK-  95.789 

LEVEL  5  RESIDUAL  NORM=  .284-010  WORK555  96.039 

LEVEL  6  RESIDUAL  NORM=  .473-009  WORK55  97.039 

LEVEL  6  RESIDUAL  NORM55  .194-009  WORK=  98.039 

LEVEL  6  RESIDUAL  NORM=  .110-009  WORK=  99.039 

******************* *enD  OF  CYCLE* ******* ************MU  =  .7787 

LEVEL  5  RESIDUAL  NORM=  .737-010  WORK=  99.289 

LEVEL  5  RESIDUAL  NORM55  .499-010  WORK=  99.539 

LEVEL  4  RESIDUAL  NORM 55  .331-010  WORK=  99.602 

LEVEL  4  RESIDUAL  NORM=  .231-010  WORK=  99.664 

LEVEL  3  RESIDUAL  NORM55  .145-010  WORK=  99.680 


Figure  3.2:  Typical  Output  for  the  PFAS  Algorithm. 

(M  =  6,  Problem  (4.1)-(4.2),  Run  #X67368) 
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rf  =  .110  10 


-9 


and 


Uf  = 


.110  10~ 
4.95 


-9-j  1/(99.039-1) 


=  .7787  . 


We  usually  only  quote  to  one  decimal  place  and  p^  to  two  decimal  places, 

since  this  is  quite  adequate  for  our  purposes. 

PFAS  computes  and  prints 


=  cll’u  ||G/|!vuinitial||Gi 


i/[wu-wuiniUall 


(3.6) 


at  the  end  of  each  cycle. 

7.  In  all  the  experiments  reported  here  the  parameters  6  and  n  (see  (2.29) 
and  (2.39))  were  given  by  6  =  .5  and  ri  =  .15.  According  to  Brandt  [77]  the  rate  of 
convergence  is  not  very  sensitive  to  changes  in  these  parameters,  and  this  was  con¬ 
firmed  in  a  few  experiments. 

In  a  few  cases,  but  never  for  6  =  .5  and  n  =  .15,  the  program  "hunted":  that 

Ml  k 

is,  the  program  went  down  from  G  to  G  ,  up  to  G  for  k  <  M,  and  then  down 
L  M 

again  to  G  instead  of  continuing  up  to  G  .  This  might  happen  several  times  before 
G  was  reached  again. 
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NUMERICAL  RESULTS  FOR  POROUS  FLOW  THROUGH  A  DAM. 


Calculations  were  performed  on  the  well-known  free  boundary  problem  describing 
the  flow  of  water  through  a  porous  dam.  The  geometry  is  shown  in  Figure  4.1.  Water 
seeps  from  a  reservoir  of  height  y^  through  a  rectangular  dam  of  width  a  to  a 
reservoir  of  height  y^.  Part  of  the  dam  is  saturated  and  the  remainder  of  the  dam  is 
dry.  The  wet  and  dry  regions  are  separated  by  an  unknown  free  boundary  which  must  be 
found  as  part  of  the  solution.  For  an  introduction  to  the  problem  see  Bear  [1972], 
or  Cryer  [1976] . 


T  y 

A  '  F 


Figure  4.1  Seepage  Through  a  Simple  Rectangular  Dam 

As  shown  by  Baiocchi  [1971]  the  problem  can  be  formulated  as  follows:  Find  u  on 
the  rectangle  fi  =  ABCF  such  that 

2 

V  u  £  1 ,  on  n, 
u  >_  0,  on  n, 

u(V^u  -  1)  —  0,  on  [4.1) 
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r 


U  =  g 


2 

(yx  -  y)  /2  on  ab  , 

2 

(y2  -  y)  /2  on  CD  , 

[y^(a  -  x)  +  y2<x)]/2a,  on  BC  , 
0,  on  DFA  , 


(4.2) 


which  is  in  the  form  (1.1). 

k-1  k-1 

This  problem  was  solved  using  PFAS,  with  and  being  injections 

(equation  (2.35))  and  with  1^  ^  defined  by  linear  interpolation  (equation  (2.52)). 
-M 

The  initial  values  of  u  were  obtained  by  interpolating  the  boundary  values  of  u 
linearly  in  the  x  direction.  A  listing  of  the  program  is  given  in  Appendix  A. 

We  considered  the  well-known  case,  y^  =  24,  y^  =  4,  and  a  =  16.  In  all 
computations  G^"  was  a  (2  +  1)  x  (3  +  1)  grid  with  =  8.  The  finest  grid  used 
was  with  (128  +1)  x  (192  +  1)  =  24897  grid  points. 

2 

To  give  the  reader  an  idea  of  the  solution,  the  solution  U  of  (2.2)  is  given 
to  four  decimal  places  in  Table  4.1. 


0 

4 

8 

12 

16 

24 

0 

0 

0 

0 

0 

20 

8 

2.5371 

0 

0 

0 

16 

32 

18. 1486 

6.7841 

0 

0 

12 

72 

47.2732 

24.9879 

7.9120 

0 

8 

128 

89.9564 

53.9823 

22.6601 

0 

4 

200 

146.5702 

94. 3247 

44. 7462 

0 

0 

— 

288 

218.0000 

148.0000 

78.0000 

8 

Table  4.1. 


2 

U  for  the  Dam  Problem 
(Run  #X34654) 


M 

The  numerical  results,  for  different  values  of  M,  and  f  =TOL=0,  are  given  in 
Table  4.2.  The  most  important  conclusions  are  that  convergence  always  occurred  and  that 

O 

the  convergence  factor  is  always  less  than  .81. 
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Run  # 

X34654 

X34654 

X34654 

X34654 

X34654 

_ 

PC  3567 

M 

2 

3 

4 

5 

7 

gm 

5«7 

9  x  13 

17  x  25 

33  x  49 

129  x  193 

rf 

0* 

4 (-17) * 

1(-13) 

1  (-8) 

1  (-10) 

1  (-7) 

o 

wf 

.404 

.607 

.726 

.813 

.778 

.81 

Execution 

Time  for  100 
Work  Units 
(Seconds) 

.114 

.428 

1.04 

3.55 

** 

PS0Ropt 

.18 

.49 

.71 

.84 

.92 

.96 

Table  4.2. 

Solution 

of  the  dam  problem  using  PFAS 

•Reached  round-off  level  before  100  work  units. 

••Required  70K  workspace  so  extended  storage  facility  invoked,  and  timing  not 
compatible . 

O 

we  now  compare  the  convergence  factors  uf  in  Table  4.2  with  those  for  other 
methods  of  solving  the  LCP  (2.2). 

M 

A  popular  method  of  solving  the  LCP  (1.3)  is  G  projected  SOR  (point  SOR  with 

projection)  which  has  also  been  called  "modified  SOR"  by  Cottle. 

When  using  GM  projected  SOR  it  is  observed  experimentally  that  the  values  of  uM 

M 

settle  down  quite  quickly  into  positive  values  and  zero  values.  Thereafter  G 

M  M  M 

projected  SOR  is  equivalent  to  using  point  SOR  on  the  subset  G+  =  {x  e  G  :  U  (x)  >  0} 

M 

Thus  the  asymptotic  convergence  factor  for  G  projected  SOR  is  in  general  equal  to 

M 

the  asymptotic  convergence  factor  for  point  SOR  on  G+.  It  is  known  (Varga  [1962, 
p.  294] )  that  for  a  region  of  area  A  and  for  the  finite  difference  equations  corres¬ 
ponding  to  the  five-point  difference  approximation  to  Laplace's  equation  with  stepsize 
h,  the  convergence  factor  for  the  optimum  choice  of  overrelaxation  parameter  u>  is 
approximated  quite  well  by 

P  (h)  =  - — - - r  -  1  .  (4.3) 

1  +  3.015[h  /A]? 
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-  *-  »  .  *+rm**%r  .  %  V-  •  1  •  ..  *  ,**■*•" 


! 


I 


In  the  present  case  we  do  not  know  the  area  of  G+l  but,  as  a  rough  guide,  the  area 
M 

of  G+  is  approximately  equal  to  the  area  of  n,  which  is  about  80%  of  the  area  of 
the  rectangle  ABCF .  Therefore,  for  our  present  purposes  the  asymptotic  convergence 
factor  for  gM  projected  SOR  with  optimum  choice  of  w  may  be  taken  to  be 


SORopt 


1  +  3.015IK / C .8  x  16  x  24)] 


-  1  i 


1  +  .172  h 


-  1 


(4.4) 


and  these  values  are  given  at  the  bottom  of  Table  4.2. 

As  Table  4.2  shows,  for  large  problems,  PFAS  is  faster  than  G11  projected  SOR.  On 
7 

G  ,  for  example,  the  increase  in  speed  (measured  in  work  units)  is  in. 96/in. 81  A  5.2. 

Against  this,  two  factors  must  be  borne  in  mind:  (1)  PFAS  is  more  complicated  and 

requires  more  overhead  per  work  unit;  (2)  PFAS  requires  somewhat  more  storage.  We 

discuss  these  two  factors  below,  but  before  doing  so  we  wish  to  emphasize  that  although 

these  factors  reduce  the  advantage  in  speed  of  PFAS,  the  measured  execution  times  for 

M 

PFAS  are  much  smaller  than  those  for  G  projected  SOR  (see  Tables  5.3  and  6.3). 

1.  Overhead. 

To  obtain  an  indication  of  the  additional  overhead  required  by  PFAS,  we  compared 

u  _Q 

execution  times  for  M  =  5.  We  first  used  PFAS  with  e  =  2.10  .  This  required 

96.156  work  units  and  took  3.40  seconds.  We  then  modified  PFAS  so  that  only  the  grid 

k  =  M  was  used  and  so  that  over-relaxation  was  used  with  the  over-relaxation  parameter 

M 

u  given  by  equation  (4.4).  We  were  thus  using  G  projected  SOR  with  a  nearly 
optimum  oi.  To  reduce  ||VuM||  to  eM  =  2.10  ®  required  146  work  units  and  took 
4.82  seconds.  Since 


(3. 40/96. 156)/(4. 82/146)  =  1.07 

we  conclude  that,  in  this  application,  the  additional  overhead  required  by  PFAS  only 

M 

increases  the  computation  time  per  G  work  unit  by  about  10%. 

2.  Storage. 

As  implemented  here,  PFAS  keeps  the  solutions  and  residuals  on  all  the  grids,  and 

therefore  requires  storage  for  2 [1+4  +4  +...]=  8/3  G  grids.  In  contrast, 

M  .  M 

G  projected  SOR  requires  storage  for  only  one  G  grid. 


M 

If  storage  is  at  a  premium,  the  residuals  on  G  need  not  be  stored  and  PFAS 

M 

requires  only  5/3  times  as  much  storage  as  G  projected  SOR.  If  u  is  stored  to 

-k  -]( 

double  precision,  but  u  and  b  are  stored  to  single  precision  for  k  <  M,  only 
4/3  times  as  much  storage  is  needed.  If  F(x)  were  not  the  constant  1,  but  a  com¬ 
plicated  function,  then  either  the  function  values  or  the  residuals  would  have  to  be 
M 

stored  for  G  projected  SOR,  and  PFAS  would  require  at  most  33%  more  storage. 

Finally,  the  PFMG  algorithm  described  in  Section  6  often  need  not  store  any  data  on 
M 

the  G  grid  (see  Section  6) . 

Another  possible  algorithm  for  solving  the  LCP  (1.3)  is  the  MBSOR  (modified  block 

SOR)  algorithm  of  Cottle  and  Sacher  [1978] .  This  algorithm  is  based  upon  the  solution 

of  a  sequence  of  "one-dimensional"  LCP's  in  much  the  same  way  that  line  SOR  is  based 

upon  solving  a  sequence  of  "one-dimensional"  equations.  We  used  MBSOR  to  solve  the 

dam  problem  (4.1),  (4.2),  for  the  case  M  =  5.  The  program  was  kindly  provided  by 

Professor  Sacher.  We  tried  a  few  values  of  the  over-relaxation  parameter  w,  and 

found  that  1.8  gave  the  best  results.  With  u  =  1.8  MBSOR  required  114  iterations 
M  ••  8 

to  reduce  ||Vu  ||  to  below  2.10  and  took  13.13  seconds.  The  following  comments 
G 

arise: 

1.  In  numerical  experiments  on  the  dam  problem,  Cottle  [1974]  found  that  MBSOR 

M 

was  about  20%  faster  than  "modified  point  SOR",  that  is,  G  projected  SOR.  This  is 

consistent  with  the  fact  that,  for  equations,  the  convergence  ratio  for  line  SOR  is 

only  faster  by  a  factor  of  /2  than  point  SOR  while  there  is  more  computation  per 

M 

iteration.  This  is  also  consistent  with  the  present  results,  since  G  projected  SOR 

—8 

required  146  iterations  to  reduce  the  residual  to  2.10  while  MBSOR  required  only 
114. 

2.  The  poor  execution  time  of  MBSOR  (13.13  seconds)  compared  to  PFAS  (3.40 
seconds)  can  be  explained  in  part  by  two  factors:  (a)  MBSOR  requires  more  computation 
per  iteration  than  is  needed  by  PFAS  for  a  single  work  unit;  (b)  the  MBSOR  program  was 
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written  for  the  case  of  general  coefficients,  while  the  PFAS  program  takes  advantage 

of  the  properties  of  the  five-point  difference  operator. 

3.  It  must  also  be  borne  in  mind  that  Cottle  and  Sacher  [1978)  found  that  MBSOR 
M 

was  three  times  as  fast  as  G  projected  SOR  for  the  journal  bearing  problem  where 
the  solution  is  zero  at  a  high  percentage  of  the  gridpoints. 

We  conclude  from  Table  4.2  and  from  the  above  discussion,  that  for  the  dam  problem 
(4.1),  (4.2)  PFAS  is  faster  than  GM  projected  SOR  and  modified  block  SOR  for  M  >_  5, 


that  is,  for  grids  of  dimension  33  *49  or  greater.  Furthermore,  we  also  conclude  that 

the  values  of  and  PgQpppt  *-n  Table  4.2  provide  a  reasonably  accurate  guide  to 

the  relative  performance  of  PFAS  and  GM  projected  SOR.  We  believe  that  PFAS  will  be 
M 

faster  than  both  G  projected  SOR  and  MBSOR  for  a  wide  range  of  problems. 

4.  For  a  grid  GM  with  N  gridpoints,  both  GM-pro jected  SOR  and  modified  block 


SOR  have  computation  times  which  are  0(N  '  ).  As  Table  4.2  shows,  the  computation 


time  for  PFAS  is  0 (N) .  Therefore,  the  performance  of  PFAS  vis-a-vis  the  other  methods 


improves  as  the  grids  become  finer. 
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ALTERNATIVE  IMPLEMENTATIONS  OF  PFAS. 


I 


In  this  section  we  discuss  alternative  implementations  of  PFAS,  the  best  of  which 
achieves  substantially  improved  performance. 

The  improvement  in  PFAS  which  might  be  possible  is  suggested  by  considering  the 
asymptotic  convergence  ratio,  saY«  f or  FAS  for  Poisson's  equation.  For  FAS,  the 

error  reduction  per  G^-sweep  is  .5.  If  each  G^-sweep  is  accompanied  by,  on  average, 

k  M 

one  G  sweep  for  1  <  k  <  M  -  1,  then  the  number  of  work  units  per  G  -sweep  is 

1  +  2_2  +  2~4  +  ...  =  4/3 
3/4 

and  the  convergence  ratio  is  (.5)  '  -  .595,  as  stated  by  Brandt  [1977,  p.  351],  In  the 

M 

present  case,  as  observed  in  Section  3,  the  average  number  of  work  units  per  G  sweep 
is 

1  +  4/[3(2n  -1)1  =  13/9  , 

so  that 

^FAS  =  (-5)9/13  =  -6188  •  t5'1) 

O  a 

This  value  of  is  observed  experimen tally.  The  worst  observed  value  of  for 

the  PFAS  results  quoted  in  Section  3  was  gf  =  .81.  Thus,  FAS  (for  equations)  is 

faster  than  PFAS  (for  LCP's)  by  a  factor  of  fcn  .81/fcn  .6188  =  2.28. 

Plausible  reasons  why  PFAS  is  slower  than  FAS  include  the  following: 

-k 

Dl :  Negative  components  of  u  . 

-k  k 

The  inequality  (2.41b)  requires  that  U  be  non-negative .  In  each  G  projected 

-k  k-1 

sweep  the  step  (2.7)  ensures  that  u  is  non-negative.  Furthermore,  if  1^  is 

-k-1 

the  injection  operator  the  initial  approximation  u  defined  by  (2.47)  is  also  non- 

1^ 

negative.  However,  (2.54)  does  not  preserve  non-negativity:  in  returning  to  G  from 
k-1  -k 

G  the  initial  approximation  u  may  have  negative  components,  and  this  is  often 
observed.  Of  course,  any  negative  components  are  removed  in  the  first  subsequent  G 
projected  sweep,  but  nevertheless  the  introduction  of  negative  components  must  retard 
convergence.  O 
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» 


D2;  Large  residuals  near  the  free  boundary. 

k  -k 

At  a  point  x  e  G  where  U  (x)  =  0  the  corresponding  residual 

Rk(x)  =  Fk{x)  -  LkUk(x)  (5.2) 

must  be  non-negative  because  of  the  inequality  (2.41a)  but  need  not  be  small.  □ 

D3:  Influence  of  the  discrete  interface. 

k  2 

The  discrete  interface  T  C  r  is  the  interface  between  the  set  of  points  where 

— k  k 

U  >0  and  the  set  of  points  where  U  =  0.  r  approximates  the  continuous  inter¬ 
face,  or  free  boundary,  F  separating  the  points  where  the  solution  u(x)  is  positive 
from  the  points  where  u(x)  is  zero. 

In  special  cases  it  may  happen  that  r  =  r  for  all  k,  in  which  case  PFAS  con¬ 
verges  as  fast  as  FAS.  An  example  is  given  by  problem  (5.3),  (5.4)  below  with  R  =  2, 

Jr 

for  which  F  is  the  line  y  =  5  -  2x;  it  is  found  experimentally  that  r  =  F  for 
k  <_  6. 

In  general,  rk  and  T  differ  by  O(h^),  and  Fk  and  Tk  1  differ  by  O(h^). 

*~k  —  k”l_ 

In  particular,  it  may  happen  that  U  (x)  >  0  while  U  (x)  =  0.  Furthermore,  near 
k-1 

F  the  residuals  may  be  less  smooth  because  of  the  projection  (2.7)  and  because  of 

k  k-1 

the  irregular  shape  of  r  and  F  .  This  introduces  errors  in  the  coarse  grid 
corrections  (2.55)  thereby  slowing  the  rate  of  convergence.  Finally,  the  injec¬ 
tion  operator  (2.35)  is  not  adequate  if  the  data  to  which  it  is  applied  is  not 

smooth .  □ 

k-1 

Multigrid  algorithms  can  often  be  speeded  up  by  modifying  the  operators  I  , 
k-1  k 

,  and  1^  We  have  tried  a  number  of  modifications  of  the  corresponding  PFAS 
subroutines  which  were  intended  to  address  the  difficulties  D1  to  D3  mentioned  above. 

The  subroutine  PFAS  in  Appendix  A  was  modified  so  as  to  facilitate  experimenta¬ 
tion.  This  was  done  by  changing  the  calls  to  the  auxiliary  subroutines  so  that  input 
"switch"  parameters  determined  which  version  of  each  subroutine  was  used. 
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In  addition,  computations  were  also  made  for  the  following  problem: 


(a)  V2 u  <_  f(x,y),  in  ft  , 

(b)  u  0,  in  n  ,  (5.3) 

(c)  u  =  g,  on  , 

where  H  =  [0,3]  *  [0,2],  and  where  f  and  g  are  chosen  so  that  the  exact  solution 
is 

2 

u  =  [cos  (x  +  y)  +2]  [max  (0;2. 5  R-Rx-y}]  .  (5.4) 

2 

Here,  R  is  a  parameter  which  is  chosen  close  to  the  value  2.  Note  that  u  e  C  ((2) 
and  u  =  0  above  the  line  y  =  R{2.5  -  x).  By  changing  the  value  of  R  we  can  force 

gridpoints  to  lie  very  close  to  the  exact  free  boundary;  this  may  be  expected  to  cause 

“  Jc  Jc 

PFAS  difficulty  because  if  U  (x)  is  positive  but  very  small  for  some  x  e  G  then 

*  '  -k 

it  will  take  PFAS  a  large  number  of  iterations  to  determine  whether  U  (x)  is  zero  or 
positive. 

The  modified  version  of  PFAS  is  called  PFASMD  and  is  listed  in  Appendix  B  as  part 
of  a  program  for  solving  the  porous  flow  problem  of  Section  4  and  problems  (5.3),  (5.4) 
PFASMD  was  used  to  compute  all  the  results  in  this  section. 

Our  first  modifications  to  the  auxiliary  subroutines  of  PFAS  were  not  very  success 
ful,  but  they  were  very  instructive  and  we  briefly  sunmarize  them.  In  all  cases,  the 
results  are  for  the  dam  problem  with  M  =  5.  (All  were  with  run  #X35519). 

Ml.  PFAS  was  modified  so  as  to  enforce  nonnegativity  of  u  immediately  after 
k-1 

returning  from  G  .  This  was  done  by  defining  ^  in  (2.55)  by 

^(uk;u*t)  «  max{0,uk}  .  (5.5) 

The  new  subroutine  was  called  INTAPR. 

This  modification  converged  slightly  faster  than  PFAS  with 
M2.  The  usual  situation  in  which  the  nonnegativity  of  u 
follows . 

-k  k  k-1  k-1 

Let  u  (x)  =  0,  where  x  e  G  but  x  /  G  .  Let  y  e  G  be  a  neighbor  of 

x,  such  that  uNy)  >  0.  It  may  then  happen  that  1  (y)  <  0.  As  a  result, 

k  k— 1  — k 

(Ik  ^w  ) (x)  may  be  negative,  and  if  so  the  updated  value  of  u  (x)  will  be  negative 


Wf  =  .803.  □ 

is  violated  is  as 


To  avoid  this,  PFAS  was  modified  by  changing  the  subroutines  SUBTRC  and  PUTU  so 


that  the  operator  I. 


k-I 


became 


f  uk(y)  if  uk  (x)  >  0  for  all  eight 


Jc-l-k 
lk  u  (y)  = 


neighbors  x  of  y  in  G  , 
0,  otherwise  . 

The  new  subroutines  were  called  PUTONN  and  SUBTNN,  respectively. 


(5.6) 


Remembering  from  (2.48)  that 


-k-l  k-l  k-l-k 
U  =  W  +  I,  u  , 
k 

we  see  from  (5.6)  and  (2.41b)  that  the  restraint  W*4  *  (y)  ^  0  is  enforced  for  every 

k-l  k  -)c 

point  y  e  G  with  a  neighbor  x  e  G  such  that  u  (x)  =  0. 

This  modification  converged  slightly  more  slowly  than  PFAS,  with  =  .817.  □ 

-M 

m3.  PFAS  was  modified  so  that  if  the  current  value  of  u  (x)  was  zero,  then 
u  (x)  was  forced  to  be  zero  for  k  <  M.  This  was  done  by  changing  the  subroutine 
RELAX.  In  effect,  (2.7)  was  followed  by  a  further  operation: 

If  k  <  M  and  uM(x,c)  =  0  then  uk's  =  0  . 

3  3 

The  new  subroutine  was  called  RELXFR. 


(5.7) 


This  modification  converged  but  much  more  slowly  than  PFAS  with  Uj.  =  .887  □ 

M4.  Brandt  [1977,  p.  378]  has  found  residual  weighting  useful  when  the  coefficients 

of  the  differential  equation  are  changing  rapidly.  We,  therefore,  changed  the  sub- 
_k-l 


routine  RESCAL  so  that  S, 


became : 


4  s£  1rk(x)  =  £  p(A)rk(x  +  Ah  )  , 


(5.8) 


where  A=  ( A^ ,  A2 )  for  integers  A^.Aj  and  the  only  nonzero  p  (A)  are 
P(0,0)  =  1  , 


p(0,l)  =  p(l,0)  =  p(0,-l)  =  p(-l,0)  =  j  . 
P(l.l)  =  P  (1,-1)  =  P(-1,1)  =  p  (-1,-1)  =  \  . 


(5.9) 


The  new  subroutine  was  called  RESCAV. 


1  2 

This  modification  cycled  between  G  and  G  ,  as  did  also  the  further  modifica- 
k-1 

tion  for  which  I  was  also  defined  by  (5.8),  (5.9). 

The  nonconvergence  of  the  modification  M4  requires  explanation,  and  this  is 
provided  by 
Lemma  5.1. 

Let  be  defined  by  (2.56).  For  1  k  M  let  U  be  the  solution  of  the 
LCP  (2.41),  where  satisfies  (2.45).  Finally,  let  ^  satisfy 


(I^_1(zk_1)  =  0)  =>  (zk_1  =  0),  for  all  zk_1  e  R  k_1 


(5.10) 


Then  for  PFAS  to  converge  it  is  necessary  that 


sjpv  -  Lk5k]  >0  , 


(5.11) 


Ik-15k  >  o  , 


.{■Wfy  -  LkDk]  =  0 


(5.12) 

(5.13) 


— k  -k  k— 1 

Proof :  Vte  apply  PFAS  by  setting  u  =  U  ,  and  forming  the  LCP  (2.41)  on  G  : 

k-l-k-1  -k-1 

L  U  <  F  , 


0k'1  >  0  , 


(*) 


(Uk'1)T(Lk"1Uk_1  -  Fk_1)  =  0  . 

-k-1  -k-1  k 

Solving  this  exactly  so  that  u  =  U  ,  we  then  return  to  G  .  Since  PFAS  converges, 

— k  — k 

the  new  value  of  u  given  by  (2.55)  must  be  equal  to  U  .  That  is, 

k  k-1  k  f-k-l  k-l-k, 

Viw  -  Vi[u  -  Tk  u  1  - 0  ' 

which,  from  (5.10),  implies  that 


uk-1  =  ik_1uk  . 

k 


Substituting  into  (*)  and  noting  (2.45)  we  obtain  (5.11)  through  (5.13). 


The  following  remarks  follow  from  Lenma  5.1. 


1.  Lemma  5.1  brings  out  an  interesting  difference  between  multigrid  methods  for 

-k  k-k 

equations  and  for  inequalities.  For  equations,  F  -  L  U  =0  and  conditions  (5.11)- 

k-1  k-1 

(5.13)  are  satisfied  for  any  reasonable  choice  of  and  1^  ,  but  this  is  not 

true  for  inequalities.  O 

2.  Since  U  solves  (2.41),  inequalities  (5.11)  and  (5.12)  will  certainly  hold 

k-1  k-1 

if  and  I  map  nonnegative  vectors  into  nonnegative  vectors.  In  particular, 

k-1  k-1 

this  will  be  the  case  if  S,  and  I,  take  linear  combinations  of  values  with 

k  k 

nonnegative  weights.  □ 


3.  if  s^-1 

k 

k-1 

and  I, 
k 

are  injections. 

then  (5.13)  is 

implied  by  (2.41c). 

□ 

4.  If  Sk_1 
k 

is  defined 

by  (5.8)  and  (5, 

.9)  while  I.k  ^ 
k 

is  injection  then 

(5.13) 

does  not  hold  in  general.  This  is  because  in  general  there  will  be  points  x,y  £  G 

k—1  — k  -k  k 

such  that  x  e  G  ,  U  (x)  >  0,  U  (y)  =  0,  y  is  a  neighbor  of  x  in  G  and 

(Fk  -  LkUk) (y)  >  0.  Then 

Ik-1Uk(x)  =  Uk(x)  >  0  , 

and 

(Sk_1(Fk  -  LkUk))(x)  >  i  (Fk  -  LkUk)(y)  >  0  , 
so  that  (5.13)  does  not  hold.  This  explains  why  the  modification  M4  of  PFAS  did  not 


converge.  LJ 

We  now  describe  two  further  modifications  of  PFAS  which  were  tried: 

M5.  Bearing  Lemma  5.1  in  mind  it  is  possible  to  introduce  weighted  sums  for  which 
(5.13)  does  hold.  One  choice  uses  weighted  residuals  only  near  the  boundary: 


where 


4S 


k-1 

k 


k 

r 


(x) 


4rk(x),  if  uk(x)  =  0  or  if  uk(y)  >  0 

for  all  eight  neighbors  y  e  G  of  x  , 
r  k  -k 

l  p(A)r  (x  +  Ah  )  signum  [u  (x  +  Ah  )]  , 

A  k  * 

otherwise 


signum  m 


f  1,  if  a  >  0  , 
|  0,  if  <i  =  0  , 
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(5.14) 


I 


and  where  the  weights  o (A)  are  as  in  (5,9).  This  was  done  by  an  appropriate  change 
in  the  subroutine  RESCAL;  the  new  subroutine  was  called  RESCL1. 

M6.  As  mentioned  in  D1  and  D3  above,  if  u  (x)  =  0  then  it  may  happen  that 

'  Jc  —  k  k  k-  1  _k 

u  (x)  =  u  (x)  +  I^_jW  (x)  is  not  zero.  It  can  be  argued  that  changes  of  u  (x) 

k 

from  or  to  zero  should  only  be  done  on  G  .  We,  therefore,  modified  the  subroutine 
INTADD  so  that  in  (2.55)  <f  was  defined  by 


^(uk(x);uNx)  ) 


uk(x),  if  uk(x)  >  0  , 
0,  otherwise  . 


(5.15) 


The  new  subroutine  was  called  INTADM.  □ 

The  modifications  M5  and  M6  are  independent,  and  we  solved  (4.1),  (4.2)  with 


M  =  5  and 

terminated 

different  combinations  of  M5  and  M6.  In  each  case 

when  ||VuM||_<  2  10  The  results  are  summarized 

G  *“ 

,  the  computations 

in  Table  5.1. 

Modifications 

- 

M5 

M6 

M5  and  M6 

Work  Units 

96.15 

126.12 

42.81 

43.76 

Execution  Time 

3.40 

4.67 

1.63 

1.76 

(Seconds ) 

O 

.815 

.854 

.623 

.623 

Table  5.1: 

Solution  of 

(4.1),  (4 

.2) 

with  M  = 

5  and 

M  —8 

c  =2,10  for  modifications  5  and  6. 
(Run  #X3S026) 


The  performance  of  PFAS  is  of  course  problem  dependent.  In  Table  5.2  we 

compare  modifications  5  and  6  for  the  problem  (5.3),  (5.4).  As  in  Table  5.1  we 

iterated  until  ||Vu^||  <  2.10  ®  on  . 

G 
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Modi fi cations 

M5 

M6 

M5  and  M6 

Work  Units 

73.62  74.32 

56.96 

65.57 

Execution  Time 
(Seconds ) 

3.09  3.24 

2.58 

3.01 

0 

hf 

.731  .738 

.669 

.704 

Table  5.2:  Solution  of  (5.3),  (5.4) 

with  M  =  5, 

R  =  32/15 

and 

M  -8 

e  *  2.10  for  modifications  5 

and  6. 

(Run  #X35563) 


We  conclude  from  the  results  given  in  Tables  5.1  and  5.2  that  the  use  of  modifica¬ 
tion  6  yields  substantial  improvements. 

Finally,  in  Table  5.3  we  extend  Table  4.2  by  comparing  the  measured  execution 

times  for  the  projected  SOR  method  and  the  best  modification  of  PFAS  (* P  defined  by 
k-1 

(5.15)  and  defined  by  injection)  for  the  dam  problem  for  various  values  of  M. 

In  each  case,  the  iterations  were  continued  until  ||Vu||  <  2  10  8. 

G 


M 

gm 

2 

5x7 

3 

9  x  13 

4 

17  x  25 

5 

33  x  49 

6 

65  x  97 

M 

G  Pro j ected 

M 

G  iterations 

19 

34 

69 

146 

295 

SOR 

Execution  Time 

(seconds) 

.02 

.09 

.60 

4.88 

39.37 

PFASMD 

M 

G  work  units 

23 

30.5 

38.7 

42.8 

45.7 

(M6) 

Execution  Time 

(seconds ) 

.04 

.12 

.41 

1.64 

6.57 

Table  ( 

i.3:  Comparison 

M 

of  G  projected 

SOR  and 

PFASMD  (modification  M6) 

M 

-8 

for  the  dam 

i  problem  with  M  = 

■  5  and 

e  =  2.10  . 

(Run  #'s  X35584  and  X35564) 

As  can  be  seen  from  Table  5.3,  PFASMD  is  better  than  projected  SOR  except  for 
very  small  grids. 
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6.  PFMG  (PROJECTED  FULL  MULTIGRID  ALGORITHM) 


In  this  section  we  describe  PFMG  (Projected  Full  Multigrid  Algorithm)  which  is  a 
modification  of  the  Full  Multigrid  Algorithm  of  Brandt.  The  flowchart  for  PFMG  is 
given  in  Figure  6.1.  PFMG  has  been  implemented  as  a  FORTRAN  subroutine  for  the  case 

2 

when  S2  is  a  rectangle  in  R  ,  and  X  is  the  Laplacian  operator.  This  subroutine 
is  listed  in  Appendix  C  as  part  of  the  program  for  solving  the  porous  flow  free 
boundary  problem  of  Section  4,  and  the  problem  (5.3),  (5.4). 

PFMG  differs  from  PFASMD  in  the  following  respects. 

I_.  Instead  of  beginning  on  GM,  one  begins  on  a  coarser  grid  GLIN  and  gradually 
works  up  to  GM. 

The  computations  begin  on  the  initial  grid  G  ,  2  =  LIN,  with  an  initial  approxi- 

-  2  -  2  i  j 

mation  u  .  u  is  computed  to  the  required  accuracy  using  grids  G  through  G  as 

in  the  PFASMD  implementation  of  PFAS,  except  that,  as  will  be  discussed  below,  the 

decision  to  move  to  a  different  grid  is  based  on  slightly  different  criteria. 

-2  -2+1 
Once  u  has  been  found  to  sufficient  accuracy,  the  initial  approximation  u  is 

obtained  from 


-2+1  ,2+1-2 
u  =  u 


(6.1) 


£  +  1  Z 

where  is  an  interpolation  operator  taking  grid  functions  on  G  into  grid 

3.+1  £+1 

functions  on  G  .It  is  known  (Brandt  [1977,  p.  377])  that  J?  should  be  more 

£+1  - Z 

accurate  than  If  in  order  to  preserve  the  smoothness  of  u  . 

£+1 

In  PFMG  is  implemented  as  a  subroutine  INTRP3  which  uses  cubic  interpolation. 

(To  use  INTRP3  we  must  have  l  2  and  so  LIN  >_  2 . )  INTRP3  is  based  upon  repeated  use 
of  the  cubic  interpolation  formulas 


f(i-)  =  [-f(-l)  +  9f  (0)  +  9f  (1 )  -  f(2)]/16  , 
f(|-)  =  (f(-l)  -  5f  (0)  +  15f  (1 )  +  5f(2)]/16  . 

-M 

Repeating  this  process,  we  finally  obtain  an  initial  approximation  u  on 
Thereafter,  the  computation  proceeds  essentially  as  in  PFASMD. 


(6.2) 

(6.3) 
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INPUT  SUBROUTINES: 


WIAXL?  - — - ^  ||  Vu  ||  '  f  ?  - — - s4  IR2(k)  =  NR2?  - — - — IR1  =  NR1? 


Flow  Chart  for  PFMG 


k— 1 

IX.  u  is  used  to  estimate  the  local  truncation  error  on  G 


-k 


Suppose  that  the  difference  approximations  are  of  order  p  and  that  u  can  be 


extended  to  a  smooth  function  on  ft.  Then  on  G 


k-1 


and 


k-1  k-l-k  .  2  .-k  k-1 

A  I,  u  =  h,  ,£u  +  t 

k  k-1 


k-1  k-k  .  2  »-k  ^  -<p+2)  k-1 

A  u  =  h^Jtu  +  2  T 


k-1  -k 

where  the  local  truncation  error  t  depends  upon  the  derivatives  of  u  . 

the  unknown  iu  we  obtain 


k-1  . 


2P  -  1 


,, k-1  k-l-k  k-1  k-k, 

[A  Ik  u  -  4Sk  A  u  ]  , 


, , .„k-l  ,rk  .k-k,  ,  k-1  k-l-k,  ,  ...k-l-k., 

[{4Sk  (b  -  A  u  )}  +  {A  Ik  u  }  -  {4Sk  b  }] 


(6.4) 

(6.5) 

Eliminating 

(6.6) 

(6.7) 


2r-  1 

In  PFMG,  the  first  {  }  in  (6.7)  is  evaluated  in  subroutine  RESSW;  the  second  {  } 

is  computed  and  added  to  the  first  using  subroutines  CORSRE  and  PUTU;  the  third  {  }  is 

evaluated  in  subroutine  RESBW  (which  is  a  minor  modification  of  RESSW);  and,  finally, 
k-1 

t  is  estimated  in  subroutine  TAUCAM.  The  estimate  (6.7)  is  not  accurate  near  the 


discrete  interface,  and  so  TAUCAM  computes  t 


k-1 


where 


k-1,  , 
tz  (x)  = 


Tk_1(x),  if  Gk_1  (x)  >  0 
0  ,  if  uk-1(x)  =  0 


(6.8) 


Because  of  the  lack  of  smoothness  of  the  solution  near  the  free  boundary,  it  is 
not  entirely  clear  what  the  value  of  p  should  be.  It  is  known  (Brezzi  and  Sacchi 

[1976])  that  the  convergence  of  the  finite  difference  approximations  is  probably  only 

112  2 
0(h  )  in  the  W  '  (SI)  norm,  and  Nitsche  [1975]  has  proved  0(h  In  h)  convergence 

in  the  infinity  norm.  However,  these  are  global  error  bounds,  while  we  are  concerned 

with  the  asymptotic  behavior  of  the  local  truncation  error  t.  Except  in  a  neighbor- 

j, 

hood  of  the  discrete  interface  T  ,  p  is  clearly  equal  to  2.  Since  the  choice  of 

l 

p  may  vary  over  ft,  we  could  perhaps  set  p  =  1  near  T  ,  but  the  values  of  T 
near  T  are  not  very  accurate  and  so,  for  simplicity,  we  have  taken  p  =  2 
everywhere. 


-39- 


can  be  used  in 


accurate  elsewhere .  In  PFMG  this  is  done  in  the  subroutine  TAUCAM  when  k  «  Z  -  1  and 
the  input  parameter  ITAU  =  1. 

Of  course,  this  is  only  meaningful  when  ||t^  is  small  compared  to  ||  Vu^  ||  ^ : 

if  the  iterations  are  continued  for  a  long  time  then  convergence  will  not  occur  because 
the  conditions  of  Lemma  5.1  will  be  violated,  but  PFMG  is  never  used  in  this  way.  In 
fact,  experience  with  equalities  indicates  that  when  T-extrapolation  is  used,  the  best 
procedure  is  to  avoid  relaxation  after  returning  for  the  last  time  to  the  finest  grid. 

IV.  As  already  mentioned,  the  logic  of  PFMG  is  more  complicated  than  that  of  PFAS; 
it  is  best  understood  by  consulting  Figure  6.1  and  Appendix  C.  Several  parameters  are 
introduced  and  this  enables  one  to  control  explicitly  the  number  of  G  projected 
sweeps  at  any  level  k,  and  the  number  of  cycles  at  level  Z.  If 

NR1  =  NR2  =  NCYC  =  NCYCM  =*  -1  , 

PREC  =  PRECM  =  0  , 

RATIO  =  1  , 

then  the  logic  of  PFMG  reduces  to  that  of  PFAS. 

We  now  describe  numerical  results  obtained  using  PFMG  to  solve  the  dam  problem 

(4.1),  (4.2).  In  all  cases,  G1  is  a  (2  +  1)  *  (3  +  1)  grid  and  LIN  «  2. 

To  control  the  iterations  we  set  NR1  =  2 ,  NR2  =  3 ,  NCYC  =  1 ,  NCYCLN  =  3 ,  and 

Z  k 

NCYCM  =  10.  The  result  is  that  in  each  cycle  on  grid  G  ,  two  G  projected  sweeps 

f  1  k 

are  carried  out  for  1  <  k  <_  L  as  we  descend  from  G  to  G  ,  and  one  G  projected 

If  Z 

sweep  is  carried  out  as  we  ascend  from  G  to  G  .  For  Z  =  LIN  up  to  three  G 

cycles  are  allowed,  so  that  a  good  initial  approximation  can  be  obtained.  For 
Z  M 

LIN  <  Z  <  M  only  one  G  cycle  is  allowed,  while  up  to  10  G  cycles  are  allowed. 

This  will  be  clearer  after  consulting  Figure  6.2  which  shows  the  output  for  M  =  4. 
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.32441+000 

GNORM  - 

.47357+000 

SOLUTION 

:  L  INFINITY  NORM  - 

.28800+003 

GNORM  » 

.43262+003 

RELATIVE  ERROR:  L  INFINITY  NORM  - 

.11264-002 

GNORM  = 

.10947-002 

LEVEL 

1  •  •  1 

4 

RESIDUAL  NORM- 

.609+000 

WORK- 

2.687  IR1— 

1  IR2 (K)—  1 

LEVEL 

4 

RESIDUAL  NORM- 

.189+000 

WORK- 

3.687  IR1- 

2  IR2(K)=  2 

GREEN 

NORM  OF  TAU-Z  - 

.646+000 

K-  3 

LEVEL 

3 

RESIDUAL  NORM- 

.953-001 

WORK- 

3.937  IR1— 

1  IR2  ( K )  —  1 

LEVEL 

3 

RESIDUAL  NORM- 

.637-001 

WORK— 

4.187  IR1— 

2  IR2 (K)»  2 

GREEN 

NORM  OF  TAU-Z  * 

.114+001 

K-  2 

LEVEL 

2 

RESIDUAL  NORM- 

.351-001 

WORK- 

4.250  IR1— 

1  IR2  (K )  —  1 

LEVEL 

2 

RESIDUAL  NORM- 

.194-001 

WORK- 

4.312  IR1— 

2  IR2(K)«  2 

GREEN 

NORM  OF  TAU-Z  - 

.108+001 

K-  1 

LEVEL 

1 

RESIDUAL  NORM- 

.913-002 

WORK- 

4.328  IR1- 

1  IR2 (K)-  1 

LEVEL 

1 

RESIDUAL  NORM- 

.984-003 

WORK- 

4.344  IR1- 

2  IR2  (K)-  2 

LEVEL 

1 

RESIDUAL  NORM- 

.615-004 

WORK- 

4.359  IR1— 

3  IR2  ( K )-  3 

LEVEL 

2 

RESIDUAL  NORM- 

.942-002 

WORK- 

4.422  IR1— 

1  IR2  (K)“  3 

LEVEL 

3 

RESIDUAL  NORM- 

.369-001 

WORK- 

4.672  IR1— 

1  IR2  (K)-  3 

LEVEL 

4 

RESIDUAL  NORM- 

.134+000 

WORK- 

5.672  IR1- 

1  IR2  (K)—  3 

GREEN 

NORM  OF  TAU-Z  - 

.860+000 

K-  3 

SOLUTION 

ERROR: 

L 

INFINITY 

NORM  - 

.48932-001 

GNORM  - 

.23107+000 

SOLUTION 

L 

INFINITY 

NORM  - 

.28800+003 

GNORM  - 

.15697+004 

RELATIVE 

ERROR: 

L 

INFINITY 

NORM  - 

.16990-003 

GNORM  - 

.14721-003 

TIME  AT  ELAPSE  IS  .1350  SECONDS 


Figure  6.2:  Typical  output  for  the  PFMG  algorithm 
(M  -  6,  Dam  Problem,  Run  #X67705) 
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Before  discussing  how  the  error  was  controlled,  it  is  necessary  to  distinguish 
between  the  goals  of  PFMG  and  PFAS.  Asymptotically,  PFMG  and  PFAS  behave  the  same , 
because  once  PFMG  has  reached  level  M  it  performs  essentially  like  PFAS.  The  purpose 

_M 

of  PFMG  is  to  obtain  quickly  an  approximation  u  which  satisfies  the  stopping 
criterion  (6.11),  namely 


To  achieve  this  we  set 

PRECM  -  1,  TOL  -  0,  ETA  -  10  , 

DELTA  -  0,  PREC  -  0,  RATIO  ■  1  . 

Finally,  we  set  »1AX  -  30,  and  WMAXM  »  40,  though  these  values  were  of  course 
never  reached. 

PFMG  includes  the  option  of  computing,  Hu*  -  ull^  and  ||u*  -  u||G,  where  u  is 
the  exact  solution.  For  the  dam  problem,  it  is  possible  to  compute  u  analytically 
using  elliptic  integrals  (Cryer  (1976))  but  this  has  not  yet  been  done:  we  therefore 
took  u  to  be  the  most  accurate  approximation  known  to  us,  namely  the  approximation 
u7  computed  in  double  precision  on  a  (128  +  1)  *  (192  +  1)  grid  as  described  in 
Section  4.  For  problem  (5.3),  (5.4)  the  exact  solution  is  given  by  (5.4). 

We  first  performed  a  number  of  experiments  with  M  «  2,3,4,  and  5: 

1.  T-extrapolation  (with  p  =  2)  gave  slightly  worse  results  for  the  dam  problem 
and  problem  (5.3),  (5.4). 

2.  In  contrast  to  our  experience  with  PFAS,  the  use  of  modification  6  had  only  a 
slight  effect. 

3.  It  was  thought  that  convergence  might  be  improved  by  multiplying  the  difference 
7u*C(x)  by  h  for  points  x  near  the  free  boundary  before  computing  ||  Vu*4  (x)  ||  _.  This 
was  implemented  as  a  subroutine  RELAX1  but  was  found  to  have  negligible  effect. 

All  the  results  given  below  are  for  the  case  of  no  T-extrapolation  (ITAU  =  0)  and 
no  modification  (NINTSW  =  NRESSW  =  1). 
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The  results  for  the  dam  problem  for  different  values  of  M  are  shown  in  Table  6.1. 


M 

2 

3 

4 

5 

GM  work  Units 

3.75 

6.75 

5.67 

5.41 

Execution  Time  (seconds) 

.009 

.053 

.131 

.349 

l|uM-u7||y||u|L 

.00374 

.00112 

.006169 

.0000623 

IISm-57||(/||u||g 

.00334 

-00109 

.000147 

.0000405 

IWum||g 

.889 

.157 

.134 

.0714 

II^Hc 

2.39 

1.25 

0.86 

0.60 

Table  6.1:  Solution  of  the  dam  problem  using  PFMG. 
(Run  #X67247 ) 


M”1 

Since  we  only  have  estimates  for  t  ,  it  is  not  possible  to  obtain  rigorous  error 
bounds.  Nevertheless,  it  is  interesting  to  apply  the  error  bounds  of  Section  2. 

Let  denote  the  vector  obtained  by  evaluating  the  solution  u(x)  on  GM.  Then, 

from  (6.4),  (1.1),  (2.2),  (2.3),  (2.13),  and  (3.1), 

M-M  M  M 
AU  <  b  +  T  , 


so  that,  from  Lemma  2.1, 


(6.13) 


On  the  other  hand,  from  Lenina  2.2, 


|uM-  uM||2if  I|pmII2!I^mII2. 

M 


For  the  dam  problem,  P  is  an  upper  triangular  matrix  with  at  most  two  nonzero 
elements  per  row,  and  ||pM||2  i  2.  Thus, 

l|uM  -  uM!|2<  |-|!7uM!|2  .  (6.14) 
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Combining  these  inequalities  we  obtain 


or,  equivalently, 


|0M-uM||2<f  [||t»||2+2||vum||2] 
T4 


!0M-uM||Gif  [I|t”||g  +  2||vGm||g] 
M 


Using  (6.8)  and  (6.9) ,  we  conclude  that 


|oM  -  GM||g<  j-  [-J  lh”_1HG  +  2  II  VuM||g]  • 

M  2 


Next,  we  note  that  for  the  dam  problem 


4  °*i ' 


2  2 

a  -  (£)  ♦  *  -055  >  14/256 


hM-16  2 


Thus,  finally,  for  the  dam  problem, 


For  example,  for  M  »  5  we  obtain,  using  Table  6.1,  that 


|U5  -  u5||g/||u5||g£  ~  J  (0.60  ♦  2  ( .  071 )  ]/  (5. 9  10  ) 


=  .00036  ; 


the  observed  value  quoted  in  Table  6.1  is  .000040. 


In  Table  6.2  we  repeat  the  computations  of  Table  6.1  for  the  problem  (5.3),  (5.4). 


.<*■<*  *•»»  \  H*<- 


M 

2 

3 

4 

5 

Gm  Work  Units 

3.75 

6.75 

5.672 

5.414 

Execution 

Time 

(seconds) 

-028 

.103 

.263 

.842 

ii  —M  -Mi, 
||U  -  u  II 

ylWl 

^  GO 

.0147 

.000985 

.000266 

.0000645 

iiuM-  0m|| 

g/N-m 

Ho 

.0147 

.00127 

.000376 

.0000956 

|,?g”I!g 

10.5 

.241 

.121 

.0764 

II  M-ill 

II  Tz  'U 

4.18 

1 .62 

1.10 

.749 

Table  6.2  Solution  of  problem  (5.3),  (5.4)  using  PFMG. 
(Run  #X67243) 


The  error  estimate  (6.19)  also  holds  for  the  problem  (5.3),  (5.4),  since  we  are 
usinq  the  Laplace  operator  on  a  rectangle  with  sides  in  the  ratio  2:3.  Applying  (6.19) 
we  obtain 

|iu5  -  u5|!g/|[u5||g<  fi  |  |  (.75)  +  2(. 076)  1/(1. 2  104)  , 

=  .0021  : 

the  observed  value  quoted  in  Table  6.2  is  .0000645. 

-M 

The  behavior  of  the  global  error  u  -  u  can  be  checked  using  Tables  6.1  and  6.2. 

From  Tab l  4.1  we  have 


i-5  -7, 

u  -  u  | 


1-2  -7, 

u  -  u 


1/3 


0000623 


00374 


1/3 


1.96 


-2 

In  Table  f.2  the  error  in  u  is  "abnormally  large".  However, 


I1  ”■*  I1 

!,u  -ulu 

1/2 

.0000645' 

'!G3-ul|_ 

.000985 

1/2 


1 


1 .96 


Those  results  strongly  suggest  that  the  global  error  is  0(h  ). 
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The  behavior  of  the  local  error  t  can  also  be  checked  using  Tables  6.1  and  6.2. 


From  Table  6.1, 

l|lT^llG/llT2llG)1/3  =  [-60/2.39]173  i  1/2'66  , 

while,  from  Table  6.2, 

[  II  Tz  II G7  II  Tz  II  G> 1/3  =  14-18/-749)1/3  =  1/2-82  , 

so  that  t  =  O(h^)  with  q  e  (.66,. 82).  This  explains  why  i-extrapolation  with  p  =  2 

did  not  reduce  the  computational  effort  for  these  problems.  The  essential  difficulty 

is  of  course  that  the  irregularity  of  the  discrete  interface  makes  it  difficult  to 

obtain  accurate  estimates  for  t.  In  fact,  T-extrapolation  with  p  =  1  was  found  to 

be  slower  than  T-extrapolation  with  p  =  2. 

Finally,  in  Table  6.3  we  repeat  the  computations  of  Table  5.3  for  a  tolerance 

eM  =  .0714,  the  value  of  ||Vu5[|  in  Table  6.1.  We  are  thus  comparing  the  performance 

G 

of  PFAS  (with  modification  6),  PFMG,  and  projected  SOR  for  comparable  errors. 


Method 

PFMG 

PFASMD  (M6) 

Projected  SOR 

Work  Units 

5.41 

9.64 

56.0 

l|vuM[|G 

.0714 

.0239 

.0695 

Execution  Time  ... 

.  349 

.440 

1.94 

(seconds) 

Table  6.3: 

Solution  of  the  dam  problem  for  M  = 

5  and  fM  =  .0714 

using  PFASMD  (modification  6),  PFMG, 

and  projected  SOR. 

(Runs  #X67247  and  4X67250) 

From  Table  6.3,  we  see  that  PFMG  is  faster  than  projected  SOR  even  when  only  low 
accuracy  is  required.  PFAS  and  PFMG  require  comparable  times,  but  PFMG  gives  much 
more  information  and  is,  therefore,  preferable.  PFMG  also  uses  fewer  work  units  than 
PFAS.  This  is  significant  because  the  number  of  work  units  used  is  independent  of 


the  computer.  Furthermore,  on  the  basis  of  experience  with  many  problems,  it  can  be 


said  that  the  number  of  work  units  used  does  not  vary  greatly  with  the  problem:  for 
most  operators  L  PFMG  requires  only  5.4  work  units. 

We  conclude  this  section  with  some  remarks  on  the  implementation  of  PFMG: 

1.  From  Table  6.3  we  see  that  the  execution  time  per  work  unit  of  PFMG  is  greater 
than  the  comparable  quantity  for  PFAS  by  a  factor 


.349  .440 

5.41  '  9.64 


1.41  . 


This  additional  overhead  is  probably  due  to  the  cubic  interpolation  used  by  J  ,  and 

could  perhaps  be  reduced  by  better  programming.  When  -C  is  complicated,  the  additional 

overhead  required  by  PFMG  is  relatively  much  less  significant:  it  is  only  with  a  very 

simple  operator  like  the  5-point  Laplacian  that  the  additional  overhead  is  so  expensive. 

M 

2.  In  PFMG  one  often  need  not  have  any  storage  for  the  finest  grid  G  -  not  even 

M 

external  storage.  The  algorithm  visits  G  only  twice:  at  the  beginning  of  the  last 
cycle  and  at  the  end  of  the  last  cycle. 

At  the  beginning  of  the  cycle,  the  following  operations  are  performed:  inter- 

M  M  M-l  M-l 

polation  (J  , )  ;  two  G  projected  sweeps:  and  residual  transfer  (I„  and  S,.  ). 

M“1  M  M 

M 

All  these  operations  can  be  made  in  one  passage  over  G  ,  in  such  a  way  that  only  four 
M 

columns  of  G  are  held  in  memory  at  one  time.  Each  time  a  new  column,  say  column  i, 

is  created  (by  interpolation),  a  relaxation  can  be  made  in  column  i-1,  then  the  second 

relaxation  can  already  be  made  in  column  i-2  and  the  residuals  from  column  i-3  can 

be  transferred  back  to  the  coarse  grid.  Column  i-4  can  simultaneously  be  discarded 

M 

(i.e.,  replaced  by  column  i) .  After  this  visit  to  G  all  the  information  is  avail- 

-M-l  -M-l  M-l  M 

able  (in  F  and  u  )  to  solve  the  G  problem  to  the  truncation  level  of  G  „ 

M 

The  final  return  to  G  (which  would  require  the  storaqe  of  the  previous  values 

M  M  M-l 

of  u  )  is  made  in  order  to  obtain  the  solution  on  G  rather  than  on  G  ,  but  it 

does  not  improve  its  pointwise  accuracy.  If  one  is  only  interested  in  knowing  some 

functionals  of  the  solution,  these  can  be  calculated  without  having  the  final  solution 

M  -M-l  M-l 

on  G  .  To  approximate  a  functional  X(U) ,  for  example,  one  computes  K{u  )  +  cv,  , 

M 


where  1  =  JC(uM)  -  $  uM  ^  is  the  final  solution  on  \  and  u 

M  w 

M  M-l  M-l 

is  the  last  solution  on  G  before  switching  back  to  G  .  Clearly,  aw  can  be 

M 

M  M-l 

calculated  during  the  above-mentioned  passage  on  G  .  Note  that  a  is  a  relative 

M 

truncation  correction1',  similar  to  xj]  ^ .  It  makes  the  approximation  3f(uM  S  +  ^ 

M  M 

M 

correct  to  the  G  truncation  level.  3f  need  not  be  a  linear  functional. 
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7.  CONCLUSIONS  AND  RECOMMENDATIONS. 

1.  Multigrid  methods  can  easily  be  adapted  to  handle  linear  complementarity 
f roDlems  arising  from  free  boundary  problems  (see  Table  4.2). 

2.  Multigrid  methods  are  superior  to  projected  SOR  and  modified  block  SOR  (see 
Tables  5.1  and  6.3,  and  Section  4). 

3.  For  high  accuracy  solutions  of  the  discrete  LCP,  one  should  use  PFASMD  with 
modification  6  (see  Tables  5.1  and  5.2). 

4.  For  solutions  which  are  accurate  to  within  truncation  error  one  should  use  PFMG, 
with  no  modifications  (see  Tables  6.1,  6.2,  and  6.3). 

Finally,  we  conclude  with  some  comments  suggesting  possible  future  applications  of 
multigrid  methods  to  comp’ ementar i ty  problems: 

1.  For  equalities,  experience  has  shown  that  multigrid  methods  are  as  efficient 
for  problems  where  X  is  nonlinear  as  for  problems  where  £  is  linear. 

2.  Experience  from  equalities  indicates  that  with  similar  efficiency  (just  a  few 
more  work  units)  one  can  solve  much  more  difficult  proL’ems,  such  as  problems  in  which 
the  coefficients  of  X  vary  by  orders  of  magnitude  (e.g.,  large  variations  in  the 
diffusivity  of  the  dam) .  In  such  cases  SOR  and  other  methods  converge  very  slowly. 

See  Aleouffe  et.  al.  (to  appear). 

3.  The  truncation  error  near  a  discrete  interface  cannot  be  reduced  by  using 
higher  order  approximations  because  the  second  derivatives  are  usually  discontinuous. 

A  good  way  to  improve  the  approximation  would  be  to  use  finer  mesh  sizes  near  the 
discrete  interface.  This  can  be  combined  very  effectively  with  the  multigrid  process 
(see  Brandt  [1979,  Section  3]).  In  fact,  a  vast  improvement  is  expected  if 
T-extrapolation  is  used  together  with  local  refinements.  Fine  levels  will  then  be 
used  only  near  the  interface. 
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======  APPX-A-PFAS  ===== 


1. 

C 

************************************************************ 

2. 

C 

3. 

C 

THIS  PROGRAM  SOLVES  THE  PROBLEM  OF  POROUS  FLOW  THROUGH  A 

4. 

C 

RECTANGULAR  DAM  OF  HEIGHT  Y1  AND  WIDTH  A. 

5. 

C 

THE  RESERVOIR  TO  THE  RIGHT  OF  THE  DAM  IS  OF  HEIGHT  Y2. 

6. 

C 

7. 

C 

WRITTEN  BY  ACHI  BRANDT  AND  COLIN  CRYER  AUGUST  1980 

8. 

c 

9. 

c 

ADDITIONAL  PARAMETERS  USED  ARE: 

10. 

c 

NX0 

THE  NUMBER  OF  GRID  INTERVALS  IN  THE  X-DIRECTION  IN 

11. 

c 

THE  COARSEST  GRID,  GRID  1. 

12. 

c 

NYO 

THE  NUMBER  OF  GRID  INTERVALS  IN  THE  Y-DIRECTION  IN 

13. 

c 

THE  COARSEST  GRID,  GRID  1. 

14. 

c 

HO 

THE  GRID  SIZE  IN  THE  COARSEST  GRID,  GRID  1 . 

15. 

c 

M 

THE  NUMBER  OF  GRIDS  TO  BE  USED. 

16. 

c 

TOL 

THE  TOLERANCE.  COMPUTATION  TERMINATES  IF  THE  RESIDUAL 

17. 

c 

ON  THE  FINEST  GRID  IS  LESS  THAN  TOL. 

18. 

c 

WMAX 

THE  MAXIMUM  NUMBER  OF  WORK  UNITS  PERMITTED  ON  THE 

19. 

c 

FINEST  GRID.  COMPUTATION  TERMINATES  WHEN  WMAX  IS  EXCEEDED 

ro 

o 

• 

c 

IN  PRACTICAL  CASES,  ONE  SETS  WMAX=30.  IN  THE  PRESENT  WORK 

21. 

c 

WE  OFTEN  SET  WMAX=100  SO  AS  TO  OBSERVE  THE  ASYMPTOTIC 

22. 

c 

BEHAVIOR  OF  THE  ALGORITHM. 

23. 

c 

MPRINT 

THE  GRID  TO  BE  PRINTED  AT  THE  END  OF  THE  COMPUTATION. 

24. 

c 

THAT  IS,  WE  PRINT  THE  MPRINT  SUBSET  OF  THE  FINAL  ANSWER 

25. 

c 

ON  THE  GRID  M. 

26. 

c 

NQSIZE 

SIZE  OF  ARRAY  Q 

27. 

c 

MUST  BE  CHANGED  FOR  LARGE  PROBLEMS  BY  EDITING  PROGRAM 

28. 

c 

=18000  FOR  DAM  PROBLEM  M=2,3,4,5,6 

29. 

c 

=70000  FOR  DAM  PROBLEM  M=7 

30. 

c 

31. 

c 

• 

CM 

ro 

c 

ALL  THE  PARAMETERS  ARE  SET  IN  THE  PROGRAM,  BUT  THEIR  VALUES 

33. 

c 

CAN  BE 

RESET  ON  THE  NAMELIST  INPUT  CARD  WHICH  IS  READ  IN 

34. 

c 

BY  THE 

PROGRAM. 

35. 

c 

THE  NAMELIST  CARD  MUST  BE  PROVIDED  AS  INPUT. 

36. 

c 

37. 

c 

THE  PROGRAM  SETS  UP  STORAGE  FOR  THE  SOLUTIONS  AND  RIGHT 

38. 

c 

HAND  SIDES. 

39. 

c 

THE  SOLUTIONS  ARE  STORED  IN  ARRAYS  1  TO  M. 

40. 

c 

THE  RIGHT  HAND  SIDES  (  OR,  SOMETIMES  THE  RESIDUALS  ) 

41. 

c 

ARE  STORED  IN  ARRAYS  M+1  TO  2*M. 

42. 

c 

43. 

c 

THIS  PROGRAM  WAS  USED  TO  COMPUTE  THE  RESULTS  IN  FIGURE  3.2 

44. 

c 

AND  TABLES  4.1  AND  4.2  OF  THE  MRC  REPORT. 

45. 

c 

46. 

c 

47. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

48. 

EXTERNAL  G, F 

49. 

COMMON 

/PRBDAT/Y1 , Y2 , A 

50. 

COMMON 

/QDAT/NQSIZE , NQERR 

51. 

NAMELIST  /INDAT/Y 1 , Y2,A,NX0 , NYO , HO , M, TOL , WMAX , MPRINT 

52. 

NQSIZE= 

■18000 

53. 

Y1=24 

54. 

Y2=4 

55. 

A=16 

56. 

NX0=4 

57. 

NY0=6 
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60 

61. 

62. 

63. 

64. 

65.  C 

66. 

67. 

68.  C 

69. 

70. 

71  . 

72. 

73.  C 

74.  C 

75. 

76.  C 

77 .  C 

80.  C 

79. 

80. 

81  . 

82. 

83.  C 

84.  C 

85. 

86.  C 

87.  C 

88.  C 

89.  C 

90.  C 
91  . 

92. 

93. 

94. 

95. 

96. 

97. 

98. 

99.  C 

100.  C 

101. 

102.  C 

103.  C 

104.  C 

105. 

106. 

107. 

108. 

109.  C 

110.  C 

111.  C 

112.  C 

113.  C 

114. 


TOL=0 . 

WMAX=30 . 

MPRINT= 1 
READ ( 5 , INDAT ) 

WRITE ( 6 , INDAT ) 

SET  TIME  TO  ZERO 
CALL  URTIMS (0.0) 

CALL  PFAS(NX0,NY0,H0,M,TOL,WMAX,G,F) 
PRINT  ELAPSED  TIME 
T=URTIMG( 'ELAPSED  TIME' ) 

CALL  SOLPRT ( M , MPRINT ) 

STOP 

END 


DOUBLE  PRECISION  FUNCTION  F(X,Y) 

DAM  PROBLEM 

THIS  SUBROUTINE  COMPUTES  THE  RIGHT  HAND  SIDE  OF  THE 
GOVERNING  POISSON  EQUATION  DEL* DEL  U=F. 

IMPLICIT  DOUBLE  PRECISION  <A-H,0-Z) 

F=  1 . 

RETURN 

END 


DOUBLE  PRECISION  FUNCTION  G(X,Y) 

DAM  PROBLEM 

THIS  SUBROUTINE  COMPUTES  THE  BOUNDARY  DATA  AND  THE 
INITIAL  APPROXIMATION  TO  THE  SOLUTION  U. 

THE  INITIAL  APPROXIMATION  IS  OBTAINED  BY  LINEAR  INTERPOLATION 
IN  THE  X-DIRECTION  BETWEEN  THE  GIVEN  BOUNDARY  DATA. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  /PRBDAT/Y1 ,Y2,A 
G1  =  .5*(Y1 -Y ) *  *2 
G2=.5*(Y2-Y)**2 
IF (  Y.GE.Y2)  G2=0 
G=(G1*(A-X)+  G2*X)/A 
RETURN 
END 


SUB  ROUTINE  PFAS ( NX0 , NY0 , HO , M , TOL , WMAX , U 1 , F ) 

THIS  SUBROUTINE  IS  THE  MAIN  MULTIGRID  SUBROUTINE. 

IT  INITIALIZES  THE  PROBLEM,  AND  REPEATEDLY  CALLS 

THE  SUBROUTINES  RELAX , RESCAL , PUTU,CORSRE, SUBTRC, AND  INTADD. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  /QDAT/NQS IZE , NQERR 
EXTERNAL  U1,F 
DIMENSION  EPS (10) 


SET  UP  ARRAYS  1  TO  M  FOR  THE  SOLUTIONS 
AND  ARRAYS  M+1  TO  2*M  FOR  THE  RIGHT  HAND  SIDES, 
AND  CHECK  THAT  Q  ARRAY  IS  LARGE  ENOUGH 
NQERR=0 


I 


3 
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115. 

116. 

117. 

118. 

119. 

120. 
121. 
122. 

123. 

124. 

125. 

126. 

127. 

128. 

129. 

130. 

131. 

132. 

133. 

134. 

135. 

136. 

137. 

138. 

139. 

140. 

141. 

142. 

143. 

144. 

145. 

146. 

147. 

148. 

149. 

150. 

151. 

152. 

153. 

154. 

155. 

156. 

157. 

158. 

159. 

160. 
161. 
162. 

163. 

164. 

165. 

166. 

167. 

168. 

169. 

170. 

171. 


DO  1  K«1,M 

K2-2** (K-1 ) 

CALL  GRDFN(K,NX0*K2+1,NY0*K2+1,H0/K2) 

1  CALL  GRDFN(K+M,NX0*K2+1,NY0*K2+1,H0/K2) 

PRINT  10 /NQSIZE 

10  FORMAT { '  SIZE  OF  Q  ARRAY  -  110) 

IF(NQERR. EQ. 0 )GOTO  12 

PRINT  11,N2ERR 

11  FORMAT ( '  ***  ERROR  IN  GRDFN  ***  ARRAY  Q  NOT  LARGE  ENOUGH  ***' , 
*  /,'  ARRAY  Q  SIZE  SHOULD  BE  AT  LEAST  ,  110) 

STOP 

1 2  CONTINUE 
C 

C 

C  INITIALIZE 

EPS ( M ) — TOL 
K=M 
WU-0 

CALL  PUTF(M,U1,0) 

CALL  PUTF(2*M, F,2 ) 

ETA=. 5 
DELTA*. 15 
C 

C  START  OF  MAIN  LOOP  IN  WHICH  ONE  MODIFIED  GAUSS-SEIDEL 

C  SWEEP  ON  GRID  K  IS  MADE. 

C 

5  ERR=1.E30 

3  ERRP=ERR 

CALL  RELAX (K, K+M, ERR) 

IF  (WU  ,LE.  0)  ERRBEG-ERR 
WU=WU+4 , *  * ( K-M ) 

WRITE (6,4 )K, ERR, WU 

4  FORMAT (’  LEVEL ' , 12 , '  RESIDUAL  NORM*' ,  D10. 3,'  WORK-',  F7.3) 

IF(ERR.LT.EPS (K) )GOTO  2 

IF  (WU.GE.WMAX) RETURN 

IF(K.EQ. 1 .OR.ERR/ERRP.LT.  ETA )GO  TO  3 
C 

C  GO  TO  COARSER  GRID 

IF(  K.NE.M  .OR.  WU.LE.3  )  GOTO  92 
FMU=0 . 0 

IF(  ERR.GT.O  )  FMU=(ERR/ERRBEG)**( 1 .D0/(WU-1 ) ) 

PRINT  91, FMU 

91  FORMAT ( '  ',  20 ('*'),' END  OF  CYCLE ', 20 ('*'), 'MU  -  ',F8.4) 

92  CONTINUE 

CALL  RESCAL (K, K+M, K+M-1 ) 

EPS (K-1 )=DELTA*ERR 
K=K- 1 

CALL  PUTU  ( K+ 1 ,  K ) 

CAIX  CORSRE ( K, K+M ) 

GOTO  5 
C 

C  GO  TO  FINER  GRID 

2  IF  (K.EQ.M) RETURN 
CALL  SUBTRC(K+1 ,K) 

CALL  INTADD(K,K+1 ) 

K-K+1 

GOTO  5 
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END 
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172. 

173. 

174. 

175. 

176. 

177. 

178. 

179. 

180. 
181. 
182. 

183. 

184. 

185. 

186. 

187. 

188. 

189. 

190. 

191. 

192. 

193. 

194. 

195. 

196. 

197. 

198. 

199. 

200. 
201. 
202. 

203. 

204. 

205. 

206. 

207. 

208. 

209. 

210. 
211. 
212. 

213. 

214. 

215. 

216. 

217. 

218. 

219. 

220. 
221. 
222. 

223. 

224. 

225. 

226. 

227. 

228. 


C 

C 

SUBROUTINE  CORSRE ( K,  KRHS ) 

C  APPLIES  THE  DIFFERENCE  OPERATOR  ON  GRID  K 

C  TO  THE  GRID  FUNCTION  IN  ARRAY  K,  AND  ADDS  THE  RESULT  TO  THE 

C  VALUES  IN  ARRAY  KRHS. 

C  KRHS  KRHS  K  K,0 

C  B  =  R  +  A  U 

C 

C  THE  RESULT  IS  STORED  IN  ARRAY  KRHS. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  Q( 18000 ) ,  1ST (200 ) , IRHS ( 200 ) 

CALL  KEY (K, 1ST, II , JJ,H) 

CALL  KEY (KRHS, IRHS ,11, JJ,H) 

11=11-1 
J1=JJ-1 
DO  1  1=2,11 
IR=IRHS ( I ) 

IO=I ST(I) 

IM=IST (1-1 ) 

IP=IST ( 1+1 ) 

DO  1  J=2,J1 

A=-Q(IR+J)-Q(IO+J+1 )-Q ( IO+J-1 )-Q ( IM+J )-Q ( IP+J ) 

1  Q ( IR+J )=-A-4 . *Q ( IO+J ) 

RETURN 

END 

C 

C 

SUBROUTINE  GRDFN ( N , IMAX , JMAX ,  HH  ) 

C  SETS  UP  ARRAY  N. 

C  IMAX  THE  DIMENSION  IN  THE  X  DIRECTION 

C  JMAX  THE  DIMENSION  IN  THE  Y  DIRECTION 

C  HH  THE  GRID  SIZE 

C  THE  ARRAY  NST  CONTAINS  THE  STARTING  ADDRESSES  OF  THE  ARRAYS. 

C  THE  ARRAY  IMX  CONTAINS  THE  MAXIMUM  ROW  NUMBERS 

C  THE  ARRAY  JMX  CONTAINS  THE  MAXIMUM  COL  NUMBERS 

C  THE  ARRAY  H  CONTAINS  THE  GRID  SIZES. 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON/GRD/NST ( 20 ) , IMX ( 20 ) ,JMX(20) ,H(20) 

COMMON  /QDAT/NQS I ZE , NQE  RR 
DATA  IQ/1/ 

NST (N )=IQ 
IMX (N ) =IMAX 
JMX (N )=JMAX 
H(N )=HH 

IQ=IQ+ IMAX* JMAX 
IF( IQ.LE.NQSIZE+1 )  RETURN 
NQERR=IQ-1 
END 
C 
C 

SUBROUTINE  INTADD ( KC , KF ) 

C  LINEARLY  INTERPOLATES  CORRECTION  ON  COARSE  GRID  KC 

C  AND  ADDS  TO  SOLUTION  ON  GRID  KF. 

C  KF  KF  KC  KF  KF 
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c 

c 

c 
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229. 

230. 

231. 

232. 

233. 

234. 

235. 

236. 

237. 

238. 

239. 

240. 

241. 

242. 

243. 

244. 

245. 

246. 

247. 

248. 

249. 

250. 

251. 

252. 

253. 

254. 

255. 

256. 

257. 

258. 

259. 

260. 
261. 
262. 

263. 

264. 

265. 

266. 

267. 

268. 

269. 

270. 

271. 

272. 

273. 

274. 

275. 

276. 

277. 

278. 

279. 

280. 
281. 
282. 

283. 

284. 


U  =  PHI (  I  W  +  U 
KC 


IMPLICIT  DOUBLE  PRECISION  (A-H.O-2) 

COMMON  Q(18000) , ISTC(200 ) , ISTF(200 ) 

CALL  KEY ( KC , ISTC , IIC , JJC , HC ) 

CALL  KEY(KF,ISTF,IIF, JJF,HF) 

DO  1  IC=2 , IIC 

IF=2*IC-1 

JF=1 

IFO=ISTF( IF) 

IFM=ISTF( IF-1  ) 

ICO=ISTC(IC) 

ICM=I STC ( IC- 1 ) 

DO  1  JC=2 , JJC 
JF=JF+2 

A=. 5* (Q ( ICO+JC )+Q ( ICO+JC- 1 ) ) 

AM=.5* (Q ( ICM+JC )+Q ( ICM+JC-1 ) ) 

Q(IFO+JF)  =  Q ( IFO+JF)+Q ( ICO+JC) 

Q ( IFM+JF)  =  Q( IFM+JF)+. 5* (Q ( ICO+JC )+Q  ( ICM+JC ) ) 

Q( IFO+JF-1 )=Q { IFO+JF-1 )+A 
1  Q ( IFM+JF- 1 )  =  Q ( IFM+JF- 1 )+. 5* ( A+AM) 

RETURN 

END 

C 

C 

SUBROUTINE  KEY (K, 1ST, IMAX , JMAX , HH ) 

C  RECOVERS  THE  INFORMATION  ABOUT  ARRAY  K  SET  UP  BY 

C  THE  SUBROUTINE  GRDFN. 

C  THE  VALUE  OF  THE  GRID  FUNCTION  AT  THE  POINT  (I,J) 

C  IS  ADDRESSED  AS  U(IST(J)+I). 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON/GRD/NST <20 ) , IMX ( 20 ) , JMX ( 20 ) , H ( 20 ) 

DIMENSION  1ST ( 1  ) 

IMAX=IMX(K) 

JMAX= JMX ( K ) 

IS=NST(K)-JMAX-1 
DO  1  1=1, IMAX 
IS=IS  +  JMAX 
1  1ST ( I )=IS 
HH=H(K) 

RETURN 

END 

C 

C 

SUBROUTINE  PUTF ( K , F , NH ) 

C  INSERTS  THE  VALUES  OF  THE  FUNCTION  F 

C  EVALUATED  AT  THE  POINTS  OF  GRID  K 

C  AND  MULTIPLIED  BY  GRIDSIZE**NH 

C  INTO  THE  ARRAY  K. 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  Q{ 18000) , 1ST (600) 

CALL  KEY  (K, 1ST, II , JJ , H ) 

H2=H**NH 
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286. 

DO  1  J=1,JJ 

287. 

X=(I-1 )*H 

288. 

Y= ( J-1 ) *H 

289. 

1  Q(IST(I )+J)=F(X,Y)*H2 

290. 

RETURN 

291 . 

END 

292. 

C 

293. 

C 

294. 

SUBROUTINE  PUTU(KF,KC) 

295. 

c 

THIS  SUBROUTINE  INJECTS  THE  SOLUTION  ON  THE  FINE  GRID 

296. 

c 

KF  INTO  THE  COARSE  GRID  KC. 

297. 

c 

KC,0  KC  KF 

298. 

c 

U  =  I  U 

299. 

c 

KF 

300. 

c 

301  . 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

302. 

COMMON  Q ( 1 8000 ) , IUF( 200 ) , IUC ( 200 ) 

303. 

CALL  KEY(KF,IUF,IIF,JJF,HF) 

304. 

CALL  KEY ( KC, IUC, IIC, JJC/HC ) 

305. 

DO  1  IC=1,IIC 

306. 

IF=2*IC-1 

307. 

IFO=IUF( IF) 

308. 

ICO=IUC(IC) 

309. 

JF=- 1 

310. 

DO  1  JC=1,JJC 

311. 

JF=JF+2 

312. 

Q( ICO+JC )  =  Q ( IFO+JF ) 

313. 

1  CONTINUE 

314. 

RETURN 

315. 

END 

316. 

c 

317. 

c 

318. 

SUBROUTINE  RELAX ( K , KRHS , ERR ) 

319. 

c 

CARRIES  OUT  ONE  MODIFIED  GAUSS-SEIDEL 

320. 

c 

SWEEP  ON  THE  GRID  K  WITH  RIGHT  HAND  SIDE  IN  ARRAY  KRHS 

321. 

c 

RETURNS  WITH  ERR=  G-NORM  OF  THE  DYNAMIC  RESIDUALS 

322. 

c 

323. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

324. 

COMMON  Q( 18000) , I ST (200) , IRHS ( 200 ) 

325. 

CALL  KEY(K, IST,II, JJ,H) 

326. 

CALL  KEY ( KRHS , IRHS, II,JJ,H) 

327. 

11=11-1 

328. 

J1=JJ-1 

329. 

ERR=0 . 

330. 

DO  1  1=2,11 

331 . 

IR=IRHS ( I ) 

332. 

IO=IST ( I ) 

333. 

IM=IST( 1-1 ) 

334. 

IP=IST (1+1 ) 

335. 

DO  1  J=2,J1 

336. 

A=Q ( IR+J )-Q { IO+J+1 )-Q ( IO+J- 1 )-Q(IM+J)-Q(IP+J) 

337. 

QT=-.25*A 

338. 

QN=MAX(0 .0 ,QT) 

339. 

ERR=ERR+ (QN-Q ( IO+J ) ) **2 

340. 

1  Q( IO+J )=QN 

341. 

ERR=SQRT { ERR ) /H 

342. 

RETURN 
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343. 

344. 

345. 

346. 

347. 

348. 

349. 

350. 

351. 

352. 

353. 

354. 

355. 

356. 

357. 

358. 

359. 

360. 

361 . 

362. 

363. 

364. 

365. 

366. 

367. 

368. 

369. 

370. 

371. 

372. 

373. 

374. 

375. 

376. 

377. 

378. 

379. 

380. 

381. 

382. 

383. 

384. 

385. 

386. 

387. 

388. 

389. 

390. 

391. 

392. 

393. 

394. 

395. 

396. 

397. 

398. 

399. 


END 

C 

C 

SUBROUTINE  RESCAL ( KF , KRF , KRC ) 

C  CALCULATES  THE  RESIDUAL  ON  GRID  KF  WITH  RIGHT  HAND  SIDE 

C  IN  ARRAY  KRF  ,  AND  INJECTS  INTO  ARRAY  KRC. 

C  BEFORE  INJECTION,  THE  RESIDUAL  IS  SCALED 

C  BY  MULTIPLYING  BY  THE  FACTOR  4  TO  TAKE  ACCOUNT  OF  THE 

C  FACT  THAT  THE  GRID  SIZE  ON  GRID  KF  IS  HALF  THE 

C  GRIDSIZE  ON  GRID  KC. 

C  KRC  KC  KRF  KF  KF 

C  R  =  4*S  (  B  -  A  U  ) 

C  KF 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  Q { 1 8000 ) , IUF( 200 ) , IRF( 200 ) , IRC ( 200 ) 

CALL  KEY ( KF , IUF , IIF , JJF , HF ) 

CALL  KEY (KRF, IRF, IIF, JJF,HF) 

CALL  KEY(KRC , IRC , IIC , JJC , HC ) 

IIC1=IIC-1 
JJC1=JJC-1 
DO  1  IC=2,IIC1 
ICR=IRC(IC) 

IF=2*IC-1 

JF=1 

IFR=IRF( IF ) 

IFO=IUF( IF) 

IFM=IUF(IF-1 ) 

IFP=IUF( IF+1 ) 

DO  1  JC=2 , JJC1 
JF=JF+2 

S=Q(IFO+JF+1 )+Q(IFO+JF-1 )+Q ( IFM+JF)+Q ( IFP+JF) 

1  Q ( ICR+JC )=4 .* (Q( IFR+JF)-S+4 . *Q ( IFO+JF) ) 

RETURN 

END 

C 

C 

SUBROUTINE  SOLPRT ( M , MPRINT ) 

C  PRINTS  THE  ARRAY  M  ON  THE  SUBARRAY  MPRINT. 

C 

IMPLICIT  DOUBLE  PRECISION  (A-HfO-Z) 

COMMON  Q( 18000) , 1ST (600) 

DIMENSION  QTEM(100) 

CALL  KEY  ( M, 1ST, II , JJ ,H ) 

INTERV=2 ** ( M-MPRINT ) 

DO  20  J=JJ, 1,-INTERV 
L=0 

DO  10  1=1,11, INTERV 

C  X  AND  Y  ARE  NOT  PRINTED  HERE,  BUT  ARE  COMPUTED  IN 

C  CASE  A  LATER  VERSION  NEEDS  THEM. 

X=(I-1 )*H 
Y=( J-1 )*H 
L=L+1 

CTEM(L)=Q(IST(I)+J) 

10  CONTINUE 

PRINT  *, (CTEM(LL) ,LL=1 ,L) 

20  CONTINUE 
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RETURN 

END 

C 

c 

SUBROUTINE  SUBTRC ( KF , KC ) 

C  THIS  SUBROUTINE  COMPUTES  THE  VALUE  INJECTED  FROM  GRID  KF  TO 

C  GRID  KC  AND  SUBTRACTS  IT  FROM  THE  SOLUTION  ON  GRID  KC. 

C  KC  KC  KC  KF 

C  W  =  U  -  I  U 

C  KF 

C 

IMPLICIT  DOUBLE  PRECISION  <A-H,0-Z) 

COMMON  Q(18000) ,IUF(200) , IUC(200) 

CALL  KEY (KF, IUF , IIF, JJF , HF ) 

CALL  KEY(KC,IUC,IIC,JJC,HC) 

DO  1  IC=1,IIC 
IF=2*IC-1 
IFO=IUF( IF) 

ICO=IUC(IC) 

JF=-1 

DO  1  JC=1,JJC 
JF=JF+2 

Q ( ICO+JC ) =Q ( ICO+JC ) -Q ( IFO+JF ) 

1  CONTINUE 
RETURN 
END 
C 
C 
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1.  C 

2.  C 

3.  C 

4.  C 

5.  C 

6.  C 

7.  C 

8.  C 

9.  C 

10.  C 

11  .  C 

12.  C 

13.  C 

14.  C 

15.  C 

16.  C 

17.  C 

18.  C 

19.  C 

20.  C 

21.  C 

22.  C 

23.  C 

24.  C 

25.  C 

26.  C 

27.  C 

28.  C 

29.  C 

30.  C 

31.  C 

32.  C 

33.  C 

34.  C 

35.  C 

36.  C 

37.  C 

38.  C 

39.  C 

40.  C 

41.  C 

42.  C 

43.  C 

44.  C 

45.  C 

46.  C 

47.  C 

48.  C 

49.  C 

50.  C 

51.  C 

52.  C 

53.  C 

54.  C 

55.  C 

56.  C 

57.  C 


****** *********************************************** ******** 

THIS  PROGRAM  SOLVES  THE  PROBLEM  OF  POROUS  FLOW  THROUGH  A 
RECTANGULAR  DAM  OF  HEIGHT  Y1  AND  WIDTH  A. 

THE  RESERVOIR  TO  THE  RIGHT  OF  THE  DAM  IS  OF  HEIGHT  Y2. 

WRITTEN  BY  ACHI  BRANDT  AND  COLIN  CRYER  AUGUST  1980 

THIS  PROGRAM  WAS  USED  TO  COMPUTE  THE  RESULTS  IN 
SECTION  5  AND  TABLE  6.4  OF  THE  MRC  REPORT. 

ADDITIONAL  PARAMETERS  USED  ARE: 

NX0  THE  NUMBER  OF  GRID  INTERVALS  IN  THE  X-DIRECTION  IN 

THE  COARSEST  GRID,  GRID  1. 

NY0  THE  NUMBER  OF  GRID  INTERVALS  IN  THE  Y-DIRECTION  IN 

THE  COARSEST  GRID,  GRID  1. 

HO  THE  GRID  SIZE  IN  THE  COARSEST  GRID,  GRID  1 . 

M  THE  NUMBER  OF  GRIDS  TO  BE  USED. 

TOL  THE  TOLERANCE.  COMPUTATION  TERMINATES  IF  THE  RESIDUAL 

ON  THE  FINEST  GRID  IS  LESS  THAN  TOL. 

WMAX  THE  MAXIMUM  NUMBER  OF  WORK  UNITS  PERMITTED  ON  THE 

FINEST  GRID.  COMPUTATION  TERMINATES  WHEN  WMAX  IS  EXCEEDED 
IN  PRACTICAL  CASES,  ONE  SETS  WMAX=30.  IN  THE  PRESENT  WORK 
WE  OFTEN  SET  WMAX=100  SO  AS  TO  OBSERVE  THE  ASYMPTOTIC 
BEHAVIOR  OF  THE  ALGORITHM. 

MPRINT  THE  GRID  TO  BE  PRINTED  AT  THE  END  OF  THE  COMPUTATION. 

THAT  IS,  WE  PRINT  THE  MPRINT  SUBSET  OF  THE  FINAL  ANSWER 
ON  THE  GRID  M. 

NQSIZE  SIZE  OF  ARRAY  Q 

MUST  BE  CHANGED  FOR  LARGE  PROBLEMS  BY  EDITING  PROGRAM 
=18000  FOR  DAM  PROBLEM  M=2,3,4,5,6 
=70000  FOR  DAM  PROBLEM  M=7 

SWITCHES 

NFGSW  =1  DAM  PROBLEM 

=2  PROBLEM  (5. 3), (5. 4). 


NINTSW  =1  INJECTION.  SUBROUTINE  INTADD 

=2  MODIFICATION  #6.  SUBROUTINE  INTADM 

CORRECTION  ONLY  ADDED  WHEN  U.NE.0.  SEE  (5.15). 
=3  MODIFICATION  #1.  SUBROUTINE  INTAPR 
PHI=MAX(0,U) 

NPUTSW  =1  INJECTION.  SUBROUTINES  PUTU  AND  SUBTRC 

=2  MODIFICATION  #2.  SUBROUTINES  PUTUNN  AND  SUBTNN. 
TRANSFER  0  IF  ANY  NEIGHBOR  ZERO. 

NRELSW  =1  NORMAL  RELAXATION.  SUBROUTINE  RELAX 
=2  MODIFICATION  #3.  SUBROUTINE  RELXFR 
VALUES  OF  U  CHANGED  ON  GRID 
K<M  ONLY  IF  U>0  ON  GRID  M. 

NRESSW  =1  INJECTION.  SUBROUTINE  RESCAL 

=2  MODIFICATION  #5.  SUBROUTINE  RESCL1 

USES  WEIGHTED  RESIDUALS  NEAR  BOUNDARY. 
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110. 
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1  13. 
114. 


C  RESIDUALS  WITH  U<0  SET  EQUAL  TO  ZERO 

C  =3  MODIFICATION  #4.  SUBROUTINE  RESCAV 

C  USES  WEIGHTED  RESIDUALS  . 

C 

C  ALL  THE  PARAMETERS  ARE  SET  IN  THE  PROGRAM,  BUT  THEIR  VALUES 

C  CAN  BE  RESET  ON  THE  NAMELIST  INPUT  CARD  WHICH  IS  READ  IN 

C  BY  THE  PROGRAM. 

C  THE  NAMELIST  CARD  MUST  BE  PROVIDED  AS  INPUT. 

C 

C  THE  PROGRAM  SETS  UP  STORAGE  FOR  THE  SOLUTIONS  AND  RIGHT 

C  HAND  SIDES. 

C  THE  SOLUTIONS  ARE  STORED  IN  ARRAYS  1  TO  M. 

C  THE  RIGHT  HAND  SIDES  (  OR,  SOMETIMES  THE  RESIDUALS  ) 

C  ARE  STORED  IN  ARRAYS  M+1  TO  2*M. 

C 

C  ********************************************************* 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

EXTERNAL  G,F 

COMMON  /PRBDAT/Y1 ,Y2,A,R 
COMMON  /QDAT/NQSIZE , NQERR 

COMMON  / SWDAT/NFGSW , NINTSW , NPUTSW , NRELSW , NRESSW 
NAMELIST  /TNDAT/Y1, Y2,A,R,NX0 , NY0 , HO ,M, TOL, WMAX,MPRINT 
*  , NFGSW, NINTSW , NPUTSW , N RELSW , NRESSW 
CHARACTER  ITITLE(80) 

C 

C  READ  IN  AND  PRINT  TITLE  CARDS 

C  FINISH  READING  TITLE  WHEN  LAST  CARD  IS  BLANK 

C  FINISH  RUN  WHEN  TITLE  CARD  IS  BLANK 

NC=0 

5  READ  10,(ITITLE(I),I=1,80) 

10  FORMAT (80A1 ) 

NC=NC+1 

PRINT  11 , (ITITLE(I ) ,1=1 ,80) 

11  FORMAT ( 1H  , 80A1 ) 

DO  12  1=1,80 

IF  ( ITITLE ( I ) ,NE . '  ' )GOTO  5 

12  CONTINUE 
IF(NC.EQ.I)  STOP 

C 

NQSIZE=1 8000 

NFGSW= 1 

NINTSW=1 

NPUTSW=1 

NRELSW=1 

NRESSW=1 

Y1  =24 

Y2=4 

A=16 

R=32. DO/15. DO 
NX0=4 
NY0=6 
H0=4  • 

M=3 

TOL=2 • D-8 
WMAX=30. 

MPRINT=1 
READ ( 5 , INDAT ) 
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WRITE ( 6 , INDAT ) 

C  PRINT  MODIFICATION  NUMBERS 

PRINT  100 

100  FORMAT (  '0  ***  THE  FOLLOWING  MODIFICATIONS  WERE  USED 

IF(NINTSW.EQ.2)  PRINT  106 
IF(NINTSW.EQ. 3 )  PRINT  101 
IF(NPUTSW.EQ.2 )  PRINT  102 
IF(NRELSW.EQ.2)  PRINT  103 
IF(NRESSW.EQ.2 )  PRINT  105 
IF ( NRESSW . EQ . 3 )  PRINT  104 


101 

FORMAT ( 

'O' 

/ 

'MODIFICATION 

NUMBER 

1'  ) 

102 

FORMAT ( 

'O' 

9 

'MODIFICATION 

NUMBER 

2') 

103 

FORMAT ( 

'O' 

9 

'MODIFICATION 

NUMBER 

3'  ) 

104 

FORMAT ( 

'O' 

9 

'MODIFICATION 

NUMBER 

4') 

105 

FORMAT  ( 

'O' 

9 

'MODIFICATION 

NUMBER 

5') 

106 

FORMAT ('O' 
PRINT  110 

9 

'MODIFICATION 

NUMBER 

6*  ) 

110 

FORMAT  ( 

1 

*************  1  ) 

C  SET  TIME  TO  ZERO 

CALL  URTIMS (0.0) 

CALL  PFASMD (NX0 ,NY0 ,H0 ,M,TOL , WMAX ,G, F ) 

C  PRINT  ELAPSED  TIME 

T=URTIMG ( 'ELAPSED  TIME' ) 

CALL  SOLPRT(M,MPRINT) 

STOP 

END 

C 

C 

DOUBLE  PRECISION  FUNCTION  F(X,Y) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  /PRBDAT/Y 1 , Y2 , A , R 

COMMON  / SWDAT/NFGSW , NINTSW , NPUTSW , NRELSW , NRESSW 
C  THIS  SUBROUTINE  COMPUTES  THE  RIGHT  HAND  SIDE  OF  THE 

C  GOVERNING  POISSON  EQUATION  DEL*DEL  U=F. 

GOTO(  1,2), NFGSW 
C 

C  DAM  PROBLEM 

1  CONTINUE 
F=  1  . 

RETURN 

C 

C  PROBLEM  OF  SECTION  5:  (5.3)  AND  (5.4) 

2  CONTINUE 
D=2.5*R 

A=DMAX 1 ( 0 .DO , D-R*X-Y ) 

B=X+Y 

C=2  * ( R**2+ 1 ) 

F= (C-2 . *A*A) *DCOS ( B)  +4*(R+1 ) *A*DSIN ( B) +2 *C 
RETURN 
END 
C 
C 

DOUBLE  PRECISION  FUNCTION  G(X,Y) 

C  THIS  SUBROUTINE  COMPUTES  THE  BOUNDARY  DATA  AND  THE 

C  INITIAL  APPROXIMATION  TO  THE  SOLUTION  U. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  /PRBDAT/Y 1,Y2, A, R 
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172. 

COMMON  /SWDAT/NFGSW , N INTSW , NPUTSW , NRELSW , NRESSW 

173. 

GOTO (  1,2) ,NFGSW 

174. 

C 

175. 

C 

DAM  PROBLEM 

176. 

C 

THE  INITIAL  APPROXIMATION  IS  OBTAINED  BY  LINEAR  INTERPOLATION 

177. 

C 

IN  THE  X-DIRECTION  BETWEEN  THE  GIVEN  BOUNDARY  DATA. 

178. 

1 

CONTINUE 

179. 

G1 =. 5* (Y1-Y ) **2 

180. 

G2=.5*(Y2-Y)**2 

181. 

IF(  Y.GE.Y2)  G2=0 

182. 

G= (G1 * ( A-X )+  G2*X )/A 

183. 

RETURN 

184. 

C 

185. 

c 

PROBLEM  OF  SECTION  5:  (5.3)  AND  (5.4) 

186. 

c 

INITIAL  APPROXIMATION  IS  A  PERTURBATION  OF  EXACT  SOLUTION 

187. 

2 

CONTINUE 

188. 

D=2.5*R 

189. 

A=DMAX1 ( 0 . DO , D-R*X-Y ) 

190. 

B=X+Y 

191  . 

G=A*  A* ( DCOS ( B ) +2 ) 

192. 

G=G+X* ( 3-X ) *Y* (2-Y ) *1 0 

193. 

RETURN 

194. 

END 

195. 

c 

196. 

c 

197. 

SUBROUTINE  PFASMD (NX0 ,NY0 , HO , M, TOL, WMAX, Ul ,F) 

198. 

c 

THIS  SUBROUTINE  IS  THE  MAIN  MULTIGRID  SUBROUTINE. 

199. 

c 

IT  INITIALIZES  THE  PROBLEM,  AND  REPEATEDLY  CALLS 

200. 

c 

THE  SUBROUTINES  RELAX, RESCAL,PUTU,CORSRE,SUBTRC, AND  INTADD. 

201  . 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

202. 

COMMON  /QDAT/NQSIZE , NQERR 

203. 

EXTERNAL  U1,F 

204. 

DIMENSION  EPS (10) 

205. 

c 

206. 

c 

207. 

c 

SET  UP  ARRAYS  1  TO  M  FOR  THE  SOLUTIONS 

208. 

c 

AND  ARRAYS  M+1  TO  2*M  FOR  THE  RIGHT  HAND  SIDES, 

209. 

c 

AND  CHECK  THAT  Q  ARRAY  IS  LARGE  ENOUGH 

210. 

NQERR=0 

211. 

DO  1  K=1 ,M 

212. 

K2=2**(K-1 ) 

213. 

CALL  GRDFN (K, NX0 *K2+1 ,NY0*K2+1 ,H0/K2 ) 

214. 

1 

CALL  GRDFN ( K+M , NX0  *K2  + 1 , NY0 *K2 + 1 , HO /K2 ) 

215. 

PRINT  10, NQSIZE 

216. 

10 

FORMAT ( '  SIZE  OF  Q  ARRAY  =  ',  110) 

217. 

IF(NQERR.EQ.O )GOTO  12 

218. 

PRINT  11, NQERR 

219. 

11 

FORMAT ('  ***  ERROR  IN  GRDFN  ***  ARRAY  Q  NOT  LARGE  ENOUGH  *' 

220. 

*  /,'  ARRAY  Q  SIZE  SHOULD  BE  AT  LEAST  =’,  110) 

221 . 

STOP 

222. 

12 

CONTINUE 

223. 

c 

224. 

c 

225. 

c 

INITIALIZE 

226. 

EPS ( M) =TOL 

227. 

K=M 

228. 

wu=o 
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230. 

231. 

232. 
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255. 

256. 

257. 

258. 

259. 

260. 
261. 
262. 

263. 

264. 

265. 
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278. 
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280. 
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282. 
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284. 

285. 


CALL  PUTF( M,U 1 , 0 ) 

CALL  PUTF( 2*M, F, 2 ) 

ETA= . 5 
DELTA*. 15 
C 

C  START  OF  MAIN  LOOP  IN  WHICH  ONE  GAUSS-SEIDEL  PROJECTED 

C  SWEEP  ON  GRID  K  IS  MADE. 

C 

5  ERR= 1 . E30 

3  ERRP=ERR 

CALL  RELSW(K,K+M,ERR) 

IF(WU  .LE.  0)  ERRBEG=ERR 
WU=WU+4.**(K-M) 

WRITE (6,4 )K, ERR,  WU 

4  FORMAT ('  LEVEL' ,12,'  RESIDUAL  NORM*' ,  D10. 3,'  WORK* ' ,  F7.3) 

IF(ERR.LT.EPS(K) )GOTO  2 

IF  (WU.GE.WMAX) RETURN 

IF(K.EQ. 1 .OR.ERR/ERRP.LT.  ETA) GO  TO  3 
C 

C  GO  TO  COARSER  GRID 

IF(  K.NE.M  .OR.  WU.LE.3  )  GOTO  92 
FMU=0 . 0 

IF(  ERR.GT.O  )  FMU= ( ERR/ERRBEG )** ( 1 . D0/( WU-1 ) ) 

PRINT  91 ,FMU 

91  FORMAT ( '  20 (  '  * '  )  , 'END  OF  CYCLE ', 20 ('*'),' MU  =  ’,F8.4) 

92  CONTINUE 

CALL  RESSW ( K, K+M , K+M- 1 ) 

EPS ( K- 1 ) =DELTA*ERR 
K=K-1 

CALL  PUTSW ( K+ 1  ,  K ) 

CALL  CORSRE ( K, K+M ) 

GOTO  5 
C 

C  GO  TO  FINER  GRID 

2  IF  (K.EQ.M) RETURN 
CALL  SUBSW( K+1 ,  K) 

CALL  INTSW(K,K+1 ) 

K=K+1 
GOTO  5 
END 
C 
C 

SUBROUTINE  CORSRE (K, KRHS ) 

C  APPLIES  THE  DIFFERENCE  OPERATOR  ON  GRID  K 

C  TO  THE  GRID  FUNCTION  IN  ARRAY  K,  AND  ADDS  THE  RESULT  TO  THE 

C  VALUES  IN  ARRAY  KRHS. 

C  KRHS  KRHS  K  K,0 

C  B  =  R  +  A  U 

C 

C  THE  RESULT  IS  STORED  IN  ARRAY  KRHS. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  Q( 1 8000 ) , 1ST (200), IRHS ( 200 ) 

CALL  KEY ( K, 1ST, II , JJ, H ) 

CALL  KEY (KRHS, IRHS , II , JJ , H) 

11=11-1 
J1=JJ-1 
DO  1  1=2,11 
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337. 

338. 

339. 

340. 

341. 

342. 


IR=IRHS ( I ) 

IO=IST(I) 

IM=IST (1-1 ) 

IP=IST ( 1+1 ) 

DO  1  J=2,J1 

A=-Q ( IR+J )-Q (IO+J+1 )-Q(IO-KT— 1 )-Q (IM+J )-Q ( IP+J ) 

1  Q ( IR+J )=-A-4 .  *Q( IO+J ) 

RETURN 

END 

C 

C 

SUBROUTINE  GRDFN (N , IMAX , JMAX ,  HH ) 

C  SETS  UP  ARRAY  N. 

C  IMAX  THE  DIMENSION  IN  THE  X  DIRECTION 

C  JMAX  THE  DIMENSION  IN  THE  Y  DIRECTION 

C  HH  THE  GRID  SIZE 

C  THE  ARRAY  NST  CONTAINS  THE  STARTING  ADDRESSES  OF  THE  ARRAYS 

C  THE  ARRAY  IMX  CONTAINS  THE  MAXIMUM  ROW  NUMBERS 

C  THE  ARRAY  JMX  CONTAINS  THE  MAXIMUM  COL  NUMBERS 

C  THE  ARRAY  H  CONTAINS  THE  GRID  SIZES. 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON/GRD/NST ( 20 ) , IMX ( 20 ) , JMX ( 20 ) ,H ( 20 ) 

COMMON  /QDAT/NQSIZE , NQERR 
DATA  IQ/1/ 

NST(N )=IQ 
IMX (N )=IMAX 
JMX (N )=JMAX 
H(N)=HH 

IQ=IQ+IMAX*JMAX 
IF( IQ.LE.NQSIZE+1 )  RETURN 
NQERR=IQ-1 
END 
C 
C 

SUBROUTINE  INTSW(KC,KF) 

C  INTERPOLATES  CORRECTION  ON  COARSE  GRID  KC 

C  AND  ADDS  TO  SOLUTION  ON  GRID  KF. 

C  KF  KF  KC  KF  KF 

C  U  PHI (  I  W  +  U  ;  U  ) 

C  KC 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  / SWDAT/NFGSW , NINTSW ,NPUTSW , NRELSW , NRESSW 
GOTO( 1 ,2,3) , NINTSW 
C 

1  CALL  INTADD ( KC , KF ) 

RETURN 

C 

2  CALL  INTADM ( KC , KF ) 

RETURN 

C 

3  CALL  INTAPR ( KC , KF ) 

RETURN 
END 

C 

c 
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.....  APX-B-PFASMD  ----- 
SUBROUTINE  INTADD ( KC , KF ) 

C  LINEARLY  INTERPOLATES  CORRECTION  ON  COARSE  GRID  KC 

C  AND  ADDS  TO  SOLUTION  ON  GRID  KF. 

C  KF  KF  KC  KF  KF 

C  U  PHI  (  I  W  +  U  ;  U  ) 

C  KC 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  Q( 18000) , ISTC (200) , ISTF(200) 

CALL  KEY(KC/ ISTC , IIC, JJC , HC) 

CALL  KEY(KF,ISTF,IIF,JJF,HF) 

DO  1  IC=2 / IIC 
IF=2*IC- 1 
JF=1 

IFO=ISTF{ IF) 

IFM=ISTF(IF-1  ) 

ICO-ISTC(IC) 

ICM=ISTC(IC-1 ) 

DO  1  JC=2, JJC 
JF-JF+2 

A= • 5* (Q (ICO+JC) +Q ( ICO+JC-1 ) ) 

AM=. 5* (Q  ( ICM+JC ) +Q ( ICM+JC- 1 ) ) 

Q( IFO+JF )  =  Q ( IFO+JF )+Q( ICO+JC) 

Q ( IFM+JF)  =  Q ( IFM+JF) +.5* (Q  ( ICO+JC )+Q ( ICM+JC) ) 

Q ( IFO+JF- 1 )=Q ( IFO+JF- 1 )+A 
1  Q ( IFM+JF- 1 )  =  Q ( IFM+JF- 1 )+. 5* (A+AM) 

RETURN 

END 

C 

SUBROUTINE  INTADM ( KC , KF ) 

C  MODIFICATION  #6. 

C  LINEARLY  INTERPOLATES  CORRECTION  ON  COARSE  GRID  KC 

C  AND  ADDS  TO  SOLUTION  ON  GRID  KF. 

C  CORRECTION  ONLY  ADDED  IF  SOLUTION  U  ON  FINE  GRID  IS 

C  NOT  ZERO.  SEE  (5.15). 

C  KF  KF  KC  KF 

C  U  =  I  U  +  U 

C  KC 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  Q( 18000) , ISTC (200), ISTF( 200) 

CALL  KEY (KC, ISTC, IIC, JJC, HC) 

CALL  KEY (KF, ISTF, IIF, JJF,HF) 

DO  1  IC=2 , IIC 

IF-2+IC-1 

JF=1 

IFO=ISTF( IF ) 

IFM=ISTF(IF-1  ) 

ICO=ISTC( IC ) 

ICM=ISTC( IC-1  ) 

DO  1  JC=2 , JJC 
JF=JF+2 

A=. 5* (Q( ICO+JC )+Q( ICO+JC-1 ) ) 

AM=. 5* (Q( ICM+JC )+Q( ICM+JC- 1 ) ) 

IF (Q( IFO+JF) .NE.O )Q( IFO+JF)  «  Q ( IFO+JF) +Q (ICO+JC ) 

IF(Q ( IFM+JF) .NE.O )Q ( IFM+JF)  =  £ ( IFM+JF)+. 5* (Q ( ICO+JC )+Q ( ICM+JC) ) 
IF(Q ( IFO+JF-1 ) .NE.O )C(IFO+JF-1 )=Q(IFO+JF-1 )+A 
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IF(Q(IFM+JF-1 ) .NE.O )Q( IFM+JF- 1 )  -  Q(IFM+JF-1 )+.5*(A+AM) 

1  CONTINUE 

RETURN 
END 
C 
C 
C 

SUBROUTINE  INTAPR(KC , KF ) 

C  MODIFICATION  #1,  PHI=MAX(0,U) 

C  LINEARLY  INTERPOLATES  CORRECTION  ON  COARSE  GRID  KC 

C  AND  ADDS  TO  SOLUTION  ON  GRID  KF. 

C  KF  KF  KC  KF  KF 

C  U  PHI (  I  W  +  U  f  U  ) 

C  KC 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  Q(18000) ,ISTC(200) ,ISTF(200) 

CALL  KEY (KC, ISTC, IIC/ JJC ,  HC) 

CALL  KEY (KF, ISTF, IIF, JJF ,HF) 

DO  1  IC=2 , IIC 

IF=2*IC-1 

JF=1 

IFO=ISTF( IF) 

IFM=ISTF( IF-1  ) 

ICO=ISTC(IC) 

ICM=ISTC( IC-1 ) 

DO  1  JC=2, JJC 
JF=JF+2 

A= . 5 * ( Q ( ICO+JC ) +Q ( ICO+JC- 1 ) ) 

AM= . 5* (Q ( ICM+JC ) +Q ( ICM+JC-1 > > 

Q( IFO+JF )  =AMAX1 (0.OD0,  Q(IFO+JF )+Q ( ICO+JC )  ) 

Q ( IFM+JF)  =AMAX1(0.0D0,  Q( IFM+JF)+. 5* (Q ( ICO+JC )+Q ( ICM+JC ) )  ) 
Q ( IFO+JF- 1 )=AMAX1 (0 .0D0 ,Q(IFO+JF-1 )+A  ) 

1  Q( IFM+JF- 1 )  =AMAX1 (O.ODO,  Q ( IFM+JF-1 ) +.5* ( A+AM)  ) 

RETURN 

END 

C 

SUBROUTINE  KEY(K, 1ST, IMAX, JMAX , HH ) 

C  RECOVERS  THE  INFORMATION  ABOUT  ARRAY  K  SET  UP  BY 

C  THE  SUBROUTINE  GRDFN. 

C  THE  VALUE  OF  THE  GRID  FUNCTION  AT  THE  POINT  (I,J) 

C  IS  ADDRESSED  AS  U(IST(J)+I). 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON/GRD/NST ( 20 ) , IMX ( 2 0 ) , JMX ( 20 ) , H { 20 ) 

DIMENSION  IST( 1 ) 

IMAX=IMX (K) 

JMAX=JMX ( K ) 

IS=NST(K)-JMAX-1 
DO  1  1=1, IMAX 
IS=IS  +  JMAX 
1  IST{ I )=IS 
HH=H(K) 

RETURN 

END 

C 

C 
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509. 

510. 
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SUBROUTINE  PUTF(K,F,NH) 

C  INSERTS  THE  VALUES  OF  THE  FUNCTION  F 

C  EVALUATED  AT  THE  POINTS  OF  GRID  K 

C  AND  MULTIPLIED  BY  GRIDSIZE**NH 

C  INTO  THE  ARRAY  K. 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  Q ( 18000 ) , 1ST (600 ) 

CALL  KEY  (K, 1ST, II, JJ , H) 

H2=H**NH 
DO  1  1=1,11 
DO  1  J=1 , JJ 
X=(I-1 )*H 
Y=(J-1 )*H 

1  Q(IST(I )+J )=F(X,Y)*H2 
RETURN 
END 
C 
C 

SUBROUTINE  PUTSW(KF,KC) 

C  THIS  SUBROUTINE  TRANSFERS  THE  SOLUTION  ON  THE  FINE  GRID 

C  KF  INTO  THE  COARSE  GRID  KC. 

C  KC,  0  KC  KF 

C  U  =  I  u 

C  KF 

C 

COMMON  / SWDAT/NFGSW , NINT SW , NPUTSW , NRELSW , NRES SW 
GOTO( 1 ,2) , NPUTSW 

1  CALL  PUTU ( KF , KC ) 

RETURN 

2  CALL  PUTUNN ( KF ,  KC ) 

RETURN 

END 

C 

C 

SUBROUTINE  PUTU(KF,KC) 

C  THIS  SUBROUTINE  INJECTS  THE  SOLUTION  ON  THE  FINE  GRID 

C  KF  INTO  THE  COARSE  GRID  KC. 

C  KC,0  KC  KF 

C  U  =  I  U 

C  KF 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  Q( 18000 ) , IUF(200) , IUC(200 ) 

CALL  KEY ( KF , IUF , IIF , JJF , HF ) 

CALL  KEY ( KC, IUC, IIC, JJC, HC ) 

DO  1  IC=1,IIC 
IF=2*IC-1 
IFO=IUF( IF) 

ICO=IUC(IC) 

JF=- 1 

DO  1  JC=1,JJC 
JF=JF+2 

Q  ( ICO+JC ) =  Q( IFO+JF) 

1  CONTINUE 
RETURN 
END 
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C 

C 

SUBROUTINE  PUTUNN ( KF , KC ) 

C  MODIFICATION  #2.  TRANSFER  0  IF  ANY  NEIGHBOR  ZERO. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  Q( 18000), IUF(200),IUC(200) 

CALL  KEY ( KF , IUF ,  IIF , JJF ,  HF ) 

CALL  KEY (KC, IUC, IIC, JJC,HC ) 

DO  1  IC=1,IIC 
IF=2*IC-1 
IFO=IUF( IF) 

ICO=IUC(IC) 

JF=- 1 

DO  1  JC=1,JJC 
JF=JF+2 

QTEMP=  Q ( IFO+JF ) 

IF  (IC.EQ.1  .OR.  IC.EQ.IIC)  GO  TO  1 
IF  (JC.EQ.1  .OR.  JC.EQ.JJC)  GO  TO  1 
IFP=IUF( IF+1 ) 

IFM=IUF( IF-1 ) 

IF(Q ( IFP+JF-1 ) . LE . 0 )  QTEMP=0 
IF(Q ( IFP+JF+1 ) .LE. 0 )  QTEMP=0 
IF ( Q ( I FP+JF ) . LE . 0 )  QTEMP=0 
IF(Q ( IFM+JF-1 ) .LE. 0 )  QTEMP=0 
IF( Q (IFM+JF+1 ) .LE. 0 )  QTEMP=0 
IF(Q ( IFM+JF) .LE. 0 )  QTEMP=0 
IF ( Q ( IFO+JF- 1  )  .LE . 0 )  QTEMP=0 
IF(Q ( IFO+JF+1 ) .LE.O )  QTEMP=0 
1  Q ( ICO+JC ) =QTEMP 

RETURN 
END 
C 
C 
C 

SUBROUTINE  RELSW(K,KRHS,ERR) 

C  CARRIES  OUT  ONE  GAUSS-SEIDEL  PROJECTED 

C  SWEEP  ON  THE  GRID  K  WITH  RIGHT  HAND  SIDE  IN  ARRAY  KRHS 

C  RETURNS  WITH  EER=  G-NORM  OF  THE  DYNAMIC  RESIDUALS 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  /SWDAT/NFGSW , NINTSW , NPUTSW, NRELSW , NRESS W 
GOTO  (1,2), NRELSW 
C 

1  CALL  RELAX (K, KRHS, ERR) 

RETURN 

C 

2  CALL  RELXFR(K, KRHS, ERR) 

RETURN 

END 

C 

C 

SUBROUTINE  RELAX (K, KRHS , ERR ) 

C  NORMAL  RELAXATION 

C  CARRIES  OUT  ONE  GAUSS-SEIDEL  PROJECTED 

C  SWEEP  ON  THE  GRID  K  WITH  RIGHT  HAND  SIDE  IN  ARRAY  KRHS 

C  RETURNS  WITH  ERR=  G-NORM  OF  THE  DYNAMIC  RESIDUALS 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
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571. 

COMMON  Q( 18000 ) , 1ST (200 ) , IRHS(200 ) 

572. 

CALL  KEY(K,IST,II,JJ,H) 

573. 

CALL  KEY (KRHS, IRHS, II, JJ,H) 

574. 

11=11-1 

575. 

J1=JJ-1 

576. 

ERR=0 . 

577. 

DO  1  1=2, I 1 

578. 

IR=IRHS ( I ) 

579. 

IO=IST ( I ) 

580. 

IM=IST (1-1  ) 

581. 

IP=IST (1+1  ) 

582. 

DO  1  J=2  ,  J1 

583. 

A=Q ( IR+J )~Q{ IO+J+1 )-Q (IO+J-1 )-Q ( IM+J )-Q ( IP+J ) 

584. 

QT=-.25*A 

585. 

QN=MAX (0.0, QT ) 

586. 

ERR=ERR+(QN-Q(IO+J) )**2 

587. 

1  Q(IO+J)=QN 

588. 

ERR=SQRT(ERR)/H 

589. 

RETURN 

590. 

END 

591. 

C 

592. 

SUBROUTINE  RELXFR(K,KRHS,ERR) 

593. 

C 

"FROZEN"  RELAXATION:  MODIFICATION  #  3 

594. 

C 

CARRIES  OUT  ONE  GAUSS-SEIDEL  PROJECTED 

595. 

C 

SWEEP  ON  THE  GRID  K  WITH  RIGHT  HAND  SIDE  IN  ARRAY  KRHS. 

596. 

c 

RETURNS  WITH  ERR=  G-NORM  OF  THE  DYNAMIC  RESIDUALS 

597. 

c 

DOES  NOT  CHANGE  VALUE  OF  U  ON  GRID  K 

598. 

c 

IF  K<M  AND  U=0  ON  GRID  M 

599. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

600. 

COMMON  Q( 18000 ) , 1ST (200 ) , IRHS ( 200 ) 

601. 

DIMENSION  ISTM(100) 

602. 

c 

ASSUMES  THAT  U  AND  RHS  ARE  STORED  ON  GRIDS  SEPARATED  BY 

603. 

M=KRHS-K 

604. 

CALL  KEY ( K , 1ST ,II,JJ,H) 

605. 

CALL  KEY ( M, ISTM, IIM, JJM , HM) 

606. 

INTERV=2  ** ( M— K ) 

607. 

CALL  KEY (KRHS, IRHS, II, JJ ,H ) 

608. 

11=11-1 

609. 

J1=JJ-1 

610. 

ERR=0 . 

61 1. 

DO  1  1=2,11 

612. 

IR=IRHS ( I ) 

613. 

IO=IST ( I ) 

614. 

IZM=ISTM( 1+INTERV* ( 1-1 ) ) 

615. 

IM=IST (1-1 ) 

616. 

IP=IST (1+1 ) 

617. 

DO  1  J=2,J1 

618. 

IF(K.EQ.M)  GO  TO  10 

619. 

QM=Q ( IZM+ 1 +INTERV* ( J-1 ) ) 

620. 

IF(2M.EQ.O )  GO  TO  1 

621. 

10 

CONTINUE 

622. 

A=Q( IR+J )-C( IO+J+1 )-Q( IO+J-1 )-Q { IM+J )-Q( IP+J) 

623. 

QT=-.25*A 

624. 

QN  =MAX ( 0 . 0 , QT ) 

625. 

ERR=ERR+ (QN-Q ( IO+J ) ) **2 

626. 

Q(IO+J)=QN 

627. 

1  CONTINUE 
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629. 

630. 

631. 
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651. 
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653. 

654. 
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656. 

657. 

658. 

659. 
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663. 

664. 

665. 
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667. 

668. 

669. 

670. 

671. 

672. 

673. 

674. 

675. 

676. 

677. 

678. 

679. 

680. 
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682. 

683. 

684. 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


C 


C 


c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


ERR=SQRT ( ERR ) /H 

RETURN 

END 

SUBROUTINE  RESSW  ( KF ,  KRF ,  KRC ) 

CALCULATES  THE  RESIDUAL  ON  GRID  KF  WITH  RIGHT  HAND  SIDE 
IN  ARRAY  KRF  ,  AND  TRANSFERS  INTO  ARRAY  KRC. 

BEFORE  TRANSFER,  THE  RESIDUAL  IS  SCALED 
BY  MULTIPLYING  BY  THE  FACTOR  4  TO  TAKE  ACCOUNT  OF  THE 
FACT  THAT  THE  GRID  SIZE  ON  GRID  KF  IS  HALF  THE 
GRIDSIZE  ON  GRID  KC. 

KRC  KC  KRF  KF  KF 
R  =  4*S  (  B  -  A  U  ) 

KF 

COMMON  /SWDAT/NFGSW , NINTSW , NPUTSW , NRELSW , NRESSW 
GOTO  (1,2, 3), NRESSW 

1  CALL  RESCAL ( KF , KRF , KRC ) 

RETURN 

2  CALL  RESCL1 (KF, KRF, KRC) 

RETURN 

3  CALL  RESCAV ( KF , KRF , KRC ) 

RETURN 

END 


SUBROUTINE  RESCAL(KF, KRF, KRC ) 

CALCULATES  THE  RESIDUAL  ON  GRID  KF  WITH  RIGHT  HAND  SIDE 
IN  ARRAY  KRF  ,  AND  INJECTS  INTO  ARRAY  KRC. 

BEFORE  INJECTION,  THE  RESIDUAL  IS  SCALED 
BY  MULTIPLYING  BY  THE  FACTOR  4  TO  TAKE  ACCOUNT  OF  THE 
FACT  THAT  THE  GRID  SIZE  ON  GRID  KF  IS  HALF  THE 
GRIDSIZE  ON  GRID  KC. 

KRC  KC  KRF  KF  KF 

R  =  4*S  (  B  -  A  U  ) 

KF 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  Q(18000), IUF( 200 ) , IRF ( 200 ) , IRC ( 200 ) 

CALL  KEY ( KF , IUF , IIF, JJF , HF ) 

CALL  KEY(KRF,IRF,IIF,JJF,HF) 

CALL  KEY (KRC, IRC , IIC, JJC, HC ) 

IIC1=IIC-1 
JJC1=JJC-1 
DO  1  IC=2,IIC1 
ICR=IRC(IC) 

IF=2*IC-1 

JF=1 

IFR=IRF( IF ) 

IFO=IUF( IF) 

IFM=IUF( IF-1 ) 

IFP=IUF( IF+1 ) 

DO  1  JC=2 , JJC1 
JF=JF+2 
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733. 
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739. 
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C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


3 


S=Q(IFO+JF+1  )  +Q  ( IFO+JF- 1  )+$  ( IFM+JF)+Q  ( TFP+JF ) 
1  Q(ICR+JC )=4 •  +  (Q  ( IFR+JF)-S+4 •  *Q  ( IFO+JF ) ) 

RETURN 

END 


SUBROUTINE  RESCL1 (KF, KRF, KRC ) 

MODIFICATION  #5  UPDATED  JUNE  23  1980 
USES  WEIGHTED  RESIDUALS  NEAR  THE  BOUNDARY 
CALCULATES  THE  RESIDUAL  ON  GRID  KF  WITH  RIGHT  HAND  SIDE 
IN  ARRAY  KRF  ,  AND  INJECTS  INTO  ARRAY  KRC. 

BEFORE  INJECTION,  THE  RESIDUAL  IS  SCALED 
BY  MULTIPLYING  BY  THE  FACTOR  4  TO  TAKE  ACCOUNT  OF  THE 
FACT  THAT  THE  GRID  SIZE  ON  GRID  KF  IS  HALF  THE 
GRIDSIZE  ON  GRID  KC. 

KRC  KC  KRF  KF  KF 
R  =4*1  (  B  -  A  U  ) 

KF 


IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  Q( 18000 ) ,IUF(200) ,IRF( 200) ,IRC(200 ) 
DIMENSION  R(9 ) 

CALL  KEY (KF, IUF, IIF, JJF,HF ) 

CALL  KEY ( KRF , IRF , IIF , JJF , HF ) 

CALL  KEY(KRC,IRC,IIC,JJC,HC) 

IIC1=IIC-1 
JJC1=JJC-1 
DO  1  IC=2 , IIC 1 
ICR=IRC(IC) 

IF=2*IC-1 

JF=1 

IFR=IRF( IF) 

IFO=IUF( IF) 

IFM=IUF( IF-1 ) 

IFP=IUF( IF+1 ) 

DO  1  JC=2 , JJC1 
JF=JF+2 

IF ( Q ( IFO+JF ) . EQ • 0 ) GOTO  2 

IF(Q ( IFP+JF+1 ) .GT.O  .AND.  Q ( IFP+JF-1 ) .GT.O  .AND. 

.AND.  Q ( IFO+JF- 1 ). GT.O  .AND. 
.AND.  2(IFM+JF-1 ) .GT.O  .AND. 
.AND.  Q ( IFP+JF  ).GT.O  )GOTO  2 


Q( IFO+JF+1 ) .GT.O 
Q(IFM+JF+1 ) .GT.O 
Q( IFM+JF  ) • GT . 0 


N=0 


DO  3  11=1,3 
I=IF+I1-2 
DO  3  J1=1 , 3 
J=JF+J1-2 


N=N+1 

IR=IRF(I) 

IO=IUF ( I ) 

IM“IUF{ 1-1 ) 

IP=IUF( 1+1 ) 

S-QUO+J+1  )+Q  ( IO+J-1  )+Q  ( IM+J  )+Q  (IP+J  ) 
S=Q  ( IR+J )+4*Q ( IO+J ) -S 
IF(Q(IO+J) .EQ.O )S=0 
R(N)=S 
CONTINUE 
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742. 

g (ICR+JC )=R( 5 )+.5* (R(2 )+R{4 )+R(6 )+R(8 )+ 

743. 

1  . 5*(R( 1 )+R(3 )+R(7 )+R(9)  )) 

744. 

GOTO  1 

745. 

2 

S=Q { IFO+JF+1 )+Q ( IFO+JF-1 )+g (IFM+JF) +g (IFP+JF) 

746. 

C ( ICR+JC ) =4 . * (Q ( IFR+JF ) -S+4 . *Q ( IFO+JF ) ) 

747. 

1 

CONTINUE 

748. 

RETURN 

749. 

END 

750. 

C 

751. 

c 

752. 

SUBROUTINE  RESC A V ( KF , KRF , KRC ) 

753. 

c 

MODIFICATION  #4 

754. 

c 

AVERAGES  RESIDUALS  OVER  NEIGHBOURING  POINTS 

755. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

756. 

COMMON  Q( 18000 ) ,IUF(200 ) , IRF ( 200 ) , IRC( 200 ) 

757. 

CALL  KEY ( KF, IUF, IIF, JJF, HF) 

758. 

CALL  KEY (KRF, IRF, IIF, JJF,HF) 

759. 

CALL  KEY ( KRC, IRC, IIC, JJC, HC ) 

760. 

c 

CLEAR  COARSE  GRID 

761. 

DO  9  1=1, IIC 

762. 

ICR=IRC( I ) 

763. 

DO  9  J=1,JJC 

764. 

9  Q ( ICR+J )=0 . 

765. 

c 

766. 

IIF  1=IIF-1 

767. 

JJF1=JJF-1 

768. 

DO  100  IF=2 , IIF 1 

769. 

IC=(IF+1 )/2 

770. 

IL=IF+1 -2  *IC 

771. 

ICR=IRC(IC) 

772. 

IFR=IRF( IF ) 

773. 

IFO=IUF (IF) 

774. 

IFM=IUF( IF-1  ) 

775. 

IFP=IUF( IF+1 ) 

776. 

DO  100  JF=2 , JJF1 

777. 

S=Q( IFO+JF+1 )+g( IFO+JF-1 ) +Q ( IFM+JF )+Q( IFP+JF) 

778. 

RES= ( Q ( IFR+JF ) -S+4 . *Q ( IFO+JF ) ) 

779. 

JC=(JF+1 )/2 

780. 

JL=JF+ 1 -2  *JC 

781. 

K=2*IL+JL+1 

782. 

GO  TO  ( 1 ,2,3,4) ,K 

783. 

1 

Q ( ICR+JC ) =g ( ICR+JC ) +RES 

784. 

GO  TO  100 

785. 

2 

RES=RES/2 

786. 

Q ( ICR+JC ) =Q ( ICR+JC ) +RES 

787. 

Q ( ICR+JC+ 1 ) =Q ( ICR+JC+ 1 ) +RES 

788. 

GO  TO  100 

789. 

3 

RES=RES/2 

790. 

Q ( ICR+JC ) =g ( ICR+JC )+RES 

791. 

ICR1=IRC( IC+ 1 ) 

792. 

Q( ICRl+JC)=Q ( ICR1+JC) +RES 

793. 

GO  TO  100 

794. 

4 

RES=RES/4 

795. 

g ( icr+jc ) =g ( icr+jc ) +res 

796. 

gdCR+JC+i  )=gdCR+JC+i  )+res 

797. 

ICR1=IRC( IC+ 1 ) 

798. 

g(ICRl+JC)=g(ICRl+JC)+RES 

APX-B-PFASMD  =«« 


799. 

800. 
801. 
802. 

803. 

804. 

805. 

806. 

807. 

808. 

809. 

810. 
811. 
812. 

813. 

814. 

815. 

816. 

817. 

818. 

819. 

820. 
821. 
822. 

823. 

824. 

825. 

826. 

827. 

828. 

829. 

830. 

831. 

832. 

833. 

834. 

835. 

836. 

837. 

838. 

839. 

840. 

841. 

842. 

843. 

844. 

845. 

846. 

847. 

848. 

849. 

850. 

851. 

852. 

853. 

854. 

855. 


Q ( ICR1  + JC+1 ) =Q ( ICR1 +JC+ 1 ) +RES 
GO  TO  100 
1 00  CONTINUE 
RETURN 
END 
C 
C 

SUBROUTINE  SOLPRT ( M, MPRINT ) 

C  PRINTS  THE  ARRAY  M  ON  THE  SUBARRAY  MPRINT. 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  Q ( 18000 ), 1ST (600) 

DIMENSION  QTEM(100) 

CALL  KEY  (M, IST,II,JJ,H) 

INTERV=2 ** ( M-MPRINT ) 

DO  20  J=JJ, 1,-INTERV 
L=0 

DO  10  1=1 ,II,INTERV 

C  X  AND  Y  ARE  NOT  PRINTED  HERE,  BUT  ARE  COMPUTED  IN 

C  CASE  A  LATER  VERSION  NEEDS  THEM. 

X= (1-1 ) *H 
Y=(J-1)*H 
L=L+ 1 

QTEM(L)=Q(IST(I)+J) 

1 0  CONTINUE 

PRINT  * , ( QTEM( LL ) ,LL=1 ,L) 

20  CONTINUE 

RETURN 
END 
C 
C 

SUBROUTINE  SUBSW(KF,KC) 

C  THIS  SUBROUTINE  COMPUTES  THE  VALUE  TRANSFERRED  FROM  GRID  KF  TO 

C  GRID  KC  AND  SUBTRACTS  IT  FROM  THE  SOLUTION  ON  GRID  KC. 

C  KC  KC  KC  KF 

C  W  =  U  -  I  U 

C  KF 

C 

COMMON  /SWDAT/NFGSW , NINTSW , NPUTSW , NRELSW , NRESSW 
GOTO (1,2) ,NPUTSW 

1  CALL  SUBTRC ( KF , KC ) 

RETURN 

2  CALL  SUBTNN ( KF , KC ) 

RETURN 

END 

C 

C 

SUBROUTINE  SUBTRC ( KF , KC ) 

C  THIS  SUBROUTINE  COMPUTES  THE  VALUE  INJECTED  FROM  GRID  KF  TO 

C  GRID  KC  AND  SUBTRACTS  IT  FROM  THE  SOLUTION  ON  GRID  KC. 

C  KC  KC  KC  KF 

C  W  =  U  -  I  U 

C  KF 

C 


IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
COMMON  Q( 18000 ) , IUF(200 ) , IUC( 200 ) 
CALL  KEY ( KF , IUF, IIF , JJF, HF) 
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856. 

857. 

858. 

859. 

860. 
861. 
862. 

863. 

864. 

865. 
866  • 

867. 

868. 

869. 

870. 

871. 

872. 

873. 

874. 

875. 

876. 

877. 

878. 

879. 

880. 
881. 
882. 

883. 

884. 

885. 

886. 

887. 

888. 

889. 

890. 

891. 

892. 

893. 

894. 

895. 

896. 

897. 

898. 

899. 

900. 

901. 

902. 

903. 

904. 


CALL  K£Y(KC,IUC,IIC,JJC,HC) 

DO  1  IC=1,IIC 
IF=2*IC- 1 
IFO=IUF{ IF) 

ICO=IUC(IC) 

JF=- 1 

DO  1  JC=1,JJC 
JF=JF+2 

Q  ( ICO+JC )  =Q  ( ICO-KIC )  -Q  ( IFO+JF ) 

1  CONTINUE 
RETURN 
END 
C 
C 

SUBROUTINE  SUBTNN ( KF , KC ) 

C  MODIFICATION  #2.  TRANSFER  0  IF  ANY  NEIGHBOR  ZERO. 


c 

THIS 

SUBROUTINE 

COMPUTES 

THE  VALUE  INJECTED  FROM  GRID 

c 

GRID 

KC  AND 

SUBTRACTS  IT 

FROM  THE  SOLUTION  ON  GRID  KC 

c 

KC 

KC 

KC 

KF 

c 

W 

■  U 

I 

U 

c 

KF 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
COMMON  Q(18000) , IUF(200 ) , IUC( 200 ) 
CALL  KEY (KF, IUF , IIF ,  JJF , HF ) 

CALL  KEY (KC, IUC, IIC, JJC, HC ) 

DO  1  IC=1 , IIC 

IF=2*IC-1 

IFO=IUF(IF) 

ICO=IUC(IC) 

JF=- 1 

DO  1  JC=1,JJC 
JF=JF+2 

QTEMP=  Q{ IFO+JF) 

IF  (IC.EQ.1  .OR.  IC.EQ.IIC)  GO  TO  1 
IF  (JC.EQ.1  .OR.  JC.EQ.JJC)  GO  TO  1 
IFP=IUF( IF+1 ) 

IFM=IUF(IF-1 ) 

IF (Q (IFP+JF-1 ) . LE. 0 )  QTEMP=0 
IF(Q ( IFP+JF+1 ) .LE. 0 )  QTEMP=0 
IF(Q ( IFP+JF) .LE. 0 )  QTEMP=0 
IF(Q ( IFM+JF- 1 ) .LE. 0 )  QTEMP=0 
IF(Q ( IFM+JF+1 ) • LE. 0 )  QTEMP=0 
IF(Q ( IFM+JF) .LE.O )  QTEMP=0 
IF(Q ( IFO+JF-1 ) .LE. 0 )  QTEMP=0 
IF(g(IFO+JF+1).LE.O)  QTEMP=0 
1  Q( ICO+JC )=Q( ICO+JC )-QTEMP 

RETURN 
END 
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1. 

C 

************************************************************* 

2. 

C 

3. 

c 

THIS  PROGRAM  SOLVES  THE  PROBLEM  OF  POROUS  FLOW  THROUGH  A 

4. 

c 

RECTANGULAR  DAM  OF  HEIGHT  VI  AND  WIDTH  A. 

5. 

c 

THE  RESERVOIR  TO  THE  RIGHT  OF  THE  DAM  IS  OF  HEIGHT  Y2. 

6. 

c 

7. 

c 

WRITTEN 

BY  ACHI  BRANDT  AND  COLIN  CRYER  AUGUST  1980 

8. 

c 

9. 

c 

THIS  PROGRAM  WAS  USED  TO  COMPUTE  THE  RESULTS  IN 

10. 

c 

SECTION 

6  OF  THE  MRC  REPORT. 

11. 

c 

12. 

c 

ADDITIONAL  PARAMETERS  USED  ARE: 

13. 

c 

NX0 

THE  NUMBER  OF  GRID  INTERVALS  IN  THE  X-DIRECTION  IN 

14. 

c 

THE  COARSEST  GRID,  GRID  1. 

15. 

c 

NY0 

THE  NUMBER  OF  GRID  INTERVALS  IN  THE  Y-DIRECTION  IN 

16. 

c 

THE  COARSEST  GRID,  GRID  1. 

17. 

c 

HO 

THE  GRID  SIZE  IN  THE  COARSEST  GRID,  GRID  1. 

18. 

c 

M 

THE  NUMBER  OF  GRIDS  TO  BE  USED. 

19. 

c 

LIN 

THE  STARTING  GRID.  LIN.GE.2 

20. 

c 

TOL 

THE  TOLERANCE 

21. 

c 

RATIO 

TOLERANCE  ON  GRID  L  IS  TOLL=TOL* RATIO **L 

22. 

c 

WMAXM 

THE  MAXIMUM  NUMBER  OF  WORK  UNITS  PERMITTED  ON  THE 

23. 

c 

FINEST  GRID.  COMPUTATION  TERMINATES  WHEN  WMAXM  IS 

24. 

c 

EXCEEDED. 

25. 

c 

WMAX 

THE  MAXIMUM  NUMBER  OF  WORK  UNITS  PERMITTED  ON  THE 

26. 

c 

GRID  L<M.  COMPUTATION  ON  GRID  L  TERMINATES  WHEN  WMAX  IS 

« 

r- 

CM 

c 

EXCEEDED. 

28. 

c 

MPRINT 

THE  GRID  TO  BE  PRINTED  AT  THE  END  OF  THE  COMPUTATION. 

29. 

c 

THAT  IS,  WE  PRINT  THE  MPRINT  SUBSET  OF  THE  FINAL  ANSWER 

30. 

c 

ON  THE  GRID  M. 

31. 

c 

NQSIZE 

SIZE  OF  ARRAY  Q 

32. 

c 

MUST  BE  CHANGED  FOR  LARGE  PROBLEMS  BY  EDITING  PROGRAM 

33. 

c 

=18000  FOR  DAM  PROBLEM  M=2,3,4,5,6 

34. 

c 

NR1 

AFTER  NR1  RELAXATIONS  ON  THE  GRID  K+1  THERE  IS  A 

» 

in 

CO 

c 

TRANSFER  TO  GRID  K. 

36. 

c 

NR2 

AFTER  A  TOTAL  NUMBER  OF  NR2  RELAXATIONS  ON  GRID  K 

37. 

c 

THERE  IS  A  TRANSFER  TO  GRID  K+1 

38. 

c 

NCYC 

MAXIMUM  NUMBER  OF  CYCLES  ON  LEVEL  L,  LIN<  L<M 

39. 

c 

NCYCLN 

MAXIMUM  NUMBER  OF  CYCLES  ON  LEVEL  LIN 

40. 

c 

NCYCM 

MAXIMUM  NUMBER  OF  CYCLES  ON  LEVEL  M 

41. 

c 

ETA 

IF  ERR.GE. ETA+ERRP  GO  TO  COARSER  GRID 

42. 

c 

DELTA 

EPS ( K- 1 ) “DELTA* (ERROR  ERR  ON  GRID  K) 

43. 

c 

PREC 

EPS ( L ) =MAX ( PREC+TAU ( L- 1 ) , TOL* RATIO* *L ) 

44. 

c 

PRECM 

EPS ( M)=MAX (PRECM+TAU ( M-1 ) ,TOL*RATIO**M) 

45. 

c 

46. 

c 

WE  CAN  ALSO  DO  TAU  EXTRAPOLATION: 

47. 

c 

ITAU 

IF  ITAU=1  DO  TAU  EXTRAPOLATION 

48. 

c 

PT 

ORDER  OF  EXTRAPOLATION 

49. 

c 

50. 

c 

SWITCHES 

51. 

c 

52. 

c 

NFGSW 

USED  IN  SUBRUTINES  F,G, SOLRED 

53. 

c 

NFGSW 

=1  DAM  PROBLEM 

54. 

c 

=2  PROBLEM  (5. 3), (5. 4). 

55. 

c 

56. 

c 

57. 

c 

NINTSW 

=1  INJECTION.  SUBROUTINE  INTADD 
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58. 

59. 

60. 
61  . 
62. 

63. 

64. 

65. 

66. 

67. 

68. 

69. 

70. 

71. 

72. 

73. 

74. 

75. 

76. 

77. 

78. 

79. 
80  . 
81  . 
82. 

83. 

84. 

85. 

86. 
87. 
88 

89. 

90. 

91 . 

92. 

93. 

94. 

95. 

96. 

97. 

98. 

99. 
100. 
101. 
102. 

103. 

104. 

105. 

106. 

107. 

108. 

109. 

110. 
111. 
112 

113. 

114. 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 


c 


=2  MODIFICATION  #6.  SUBROUTINE  INTADM 

CORRECTION  ONLY  ADDED  WHEN  U.NE.0.  SEE  (5.15). 

NRESSW  =1  INJECTION.  SUBROUTINE  RESCAL 

=2  MODIFICATION  #5.  SUBROUTINE  RESCL1 

USES  WEIGHTED  RESIDUALS  NEAR  BOUNDARY. 
RESIDUALS  WITH  U<0  SET  EQUAL  TO  ZERO 


ALL  THE  PARAMETERS  ARE  SET  IN  THE  PROGRAM,  BUT  THEIR  VALUES 
CAN  BE  RESET  ON  THE  NAMELIST  INPUT  CARD  WHICH  IS  READ  IN 
BY  THE  PROGRAM. 

THE  NAMELIST  CARD  MUST  BE  PROVIDED  AS  INPUT. 

THE  PROGRAM  SETS  UP  STORAGE  FOR  THE  SOLUTIONS  AND  RIGHT 
HAND  SIDES. 

THE  SOLUTIONS  ARE  STORED  IN  ARRAYS  1  TO  M. 

THE  RIGHT  HAND  SIDES  (  OR,  SOMETIMES  THE  RESIDUALS  ) 

ARE  STORED  IN  ARRAYS  M+1  TO  2*M. 

THE  EXACT  SOLUTION  (WHEN  KNOWN)  IS  STORED  IN  GRID  NGRSOL 
THE  VALUES  OF  TAU  ARE  STORED  IN  GRIDS  2M+1  TO  3M-1 

********************************************************* 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

EXTERNAL  G, F 

COMMON  /PRBDAT/Yl , Y2 , A, R 
COMMON  /QDAT/NQSIZE, NQERR 
COMMON  /SOLTAU/M, NGRSOL, PT 
COMMON  /SWDAT/NFGSW , NINTSW , NRESSW 

NAMELIST  /INDA T/NX0,NY0, HO, M, LIN, NR1,NR2, ETA, DELTA 

*  , TOL , RATIO , PREC , PRECM , WMAX , WMAXM , NC YC , NCYCLN , NCYCM , ITAU , PT , 

*  MPRINT, Y1 ,Y2 ,A,R 

*  ,NFGSW, NINTSW, NRELSW, NRESSW 
CHARACTER  ITITLE(80) 

READ  IN  AND  PRINT  TITLE  CARDS 
FINISH  READING  TITLE  WHEN  LAST  CARD  IS  BLANK 
FINISH  RUN  WHEN  TITLE  CARD  IS  BLANK 
PRINT  18 
18  FORMAT ( 1  Hi ) 

NC  =  0 

5  READ  1 0 , ( ITITLE (I), 1=1, 80) 

10  FORMAT ( 80A1 ) 

NC=NC+1 

PRINT  11,  (ITITLE (I), 1=1, 80) 

11  FORMAT ( 1H  , 80A1 ) 

DO  12  1=1,80 

IF  ( ITITLE (I ) ,NE.'  ' )GOTO  5 

U  CONTINUE 

IF(NC.EQ.I)  STOP 

NQSIZE=18000 
NFGSW= 1 
NINTSW= 1 
NRESSW=1 
Y1=24 
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115. 

Y2=4 

116. 

A=16 

117. 

R=32. D0/15. DO 

118. 

NX0=2 

119. 

NY0=3 

120. 

H0=8 . 

121. 

M=6 

122. 

LIN=2 

123. 

NR1  =2 

124. 

NR2=3 

125. 

ETA= 10. 

126. 

DELTA=0 

127. 

TOL=0 

128. 

RATI0=1 

129. 

PREC=0 

130. 

PRECM=1 

131. 

WMAX=30. 

132. 

WMAXM=40 

133. 

NCYC=1 

134. 

NCYCLN=3 

135. 

NCYCM=1 0 

136. 

ITAU=0 

137. 

PT  =2 

138. 

MPRINT=2 

139. 

READ ( 5 , INDAT ) 

140. 

WRITE ( 6 , INDAT ) 

141. 

C 

PRINT  MODIFICATION  NUMBERS 

142. 

PRINT  100 

143. 

100 

FORMAT (  '0  ***  THE  FOLLOWING  MODIFICATIONS  WERE  USED  *** 

144. 

IF (NINTSW.EQ. 2 )  PRINT  106 

145. 

IF(NINTSW.EQ.3)  PRINT  101 

146. 

IF(NRELSW.E£.2)  PRINT  103 

147. 

IF(NRESSW.EQ.2)  PRINT  105 

148. 

IF(NRESSW.EQ. 3 )  PRINT  104 

149. 

101 

FORMAT ('O'.  ’MODIFICATION  NUMBER  1 ' ) 

150. 

103 

FORMAT ('O',  'MODIFICATION  NUMBER  3' ) 

151. 

104 

FORMAT ('O',  ' MODIFICATION  NUMBER  4 ' ) 

152. 

105 

FORMAT ('O',  'MODIFICATION  NUMBER  5') 

153. 

106 

FORMAT ('O',  ' MODIFICATION  NUMBER  6 ' ) 

154. 

PRINT  110 

155. 

110 

FORMAT (  1  *************  1 ) 

156. 

C 

SET  TIME  TO  ZERO 

157. 

CALL  URTIMS (0.0) 

158. 

CALL  PFMG (NX0 ,NY0 , HO , LIN ,NR1 ,NR2 , ETA , DELTA 

159. 

• 

, TOL , RATIO , PREC , PRECM , WMAX , WMAXM ,  NC  YC , NCYCLN , NCYCM , ITAU , 

160. 

• 

MPRINT/G, F) 

161 . 

T=URTIMG ( ' ELAPSE  * ) 

162. 

19 

FORMAT ( 1H0,  ’  GRID-M  SOLUTION',//) 

163. 

PRINT  19 

164. 

CALL  SOLPRT(M,MPRINT ) 

165. 

PRINT  20 

166. 

20 

FORMAT ( 1 H 1 ,  '  GRID-7  SOLUTION',//) 

167. 

CALL  SOLPRT (NGRSOL,MPRINT) 

168. 

STOP 

169. 

END 
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172. 

173. 

174. 

175. 

176. 

177. 

178. 

179. 

180. 
181. 
182. 

183. 

184. 

185. 

186. 

187. 

188. 

189. 

190. 
191  . 

192. 

193. 

194. 

195. 

196. 

197. 

198. 

199. 

200. 
201. 
202. 

203. 

204. 

205. 

206. 

207. 

208. 

209. 

210. 
211. 
212. 

213. 

214. 

215. 

216. 

217. 

218. 

219. 

220. 
221. 
222. 

223. 

224. 

225. 

226. 

227. 

228. 


DOUBLE  PRECISION  FUNCTION  F(X,Y) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

COMMON  /PRBDAT/Y 1,Y2,A,R 
COMMON  /SWDAT/N FGSW , NINTSW , NRESSW 
C  THIS  SUBROUTINE  COMPUTES  THE  RIGHT  HAND  SIDE  OF  THE 

C  GOVERNING  POISSON  EQUATION  DEL*DEL  U=F. 

GOTO(  1,2) , NFGSW 
C 

C  DAM  PROBLEM 

1  CONTINUE 
F=1 . 

RETURN 

C 

C  PROBLEM  OF  SECTION  5:  (5.3)  AND  (5.4) 

2  CONTINUE 
D=2 . 5*R 

A=DMAX1 (0.D0,D-R*X-Y) 

B=X+Y 

C=2* (R**2+1 ) 

F= (C-2 . *A* A) *DCOS ( B)  +4*(R+1 )*A*DSIN ( B)+2*C 
RETURN 
END 
C 
C 

DOUBLE  PRECISION  FUNCTION  G(X,Y) 

C  THIS  SUBROUTINE  COMPUTES  THE  BOUNDARY  DATA  AND  THE 

C  INITIAL  APPROXIMATION  TO  THE  SOLUTION  U. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  /PRBDAT/Y 1 ,Y2,A,R 
COMMON  /SWDAT/NFGSW, NINTSW, NRESSW 
GOTO (  1,2), NFGSW 
C 

C  DAM  PROBLEM 

C  THE  INITIAL  APPROXIMATION  IS  OBTAINED  BY  LINEAR  INTERPOLATION 

C  IN  THE  X-DIRECTION  BETWEEN  THE  GIVEN  BOUNDARY  DATA. 

1  CONTINUE 
G1=.5*(Y1-Y)**2 
G2=. 5*(Y2-Y)**2 
IF(  Y.GE.Y2 )  G2=0 
G=(G1*(A-X)+  G2*X )/A 
RETURN 

C 

C  PROBLEM  OF  SECTION  5:  (5.3)  AND  (5.4) 

C  INITIAL  APPROXIMATION  IS  A  PERTURBATION  OF  EXACT  SOLUTION 

2  CONTINUE 
D=2 • 5*R 

A=DMAX1 ( 0 . DO , D-R*X-Y ) 

B=X+Y 

G=A*  A* ( DCOS ( B ) +2 ) 

G=G+X* (3-X ) *Y* (2-Y ) *10 
RETURN 
END 
C 

SUBROUTINE  PFMG { NX0 , NY 0 , HO , LIN , NR 1 , NR2 , ETA , DELTA 

*  , TOL , RATIO , PREC , PRECM , WMAX , WMAXM , NCYC , NCYCLN , NCYCM , I TAU , 

*  MPRINT,U1 ,F) 

C  THIS  SUBROUTINE  IS  THE  MAIN  FULL  MULTIGRID  SUBROUTINE. 
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229. 

230. 

231. 

232. 

233. 

234. 

235. 

236. 

237. 

238. 

239. 

240. 

241. 

242. 

243. 

244. 

245. 

246. 

247. 

248. 

249. 

250. 

251. 

252. 

253. 

254. 

255. 

256. 

257. 

258. 

259. 

260. 
261. 
262. 

263. 

264. 

265. 

266. 

267. 

268. 

269. 

270. 

271. 

272. 

273. 

274. 

275. 

276. 

277. 

278. 

279. 

280. 
281. 
282. 

283. 

284. 

285. 


C  IT  INITIALIZES  THE  PROBLEM,  AND  REPEATEDLY  CALLS 

C  THE  SUBROUTINES  RELAX, RESCAL, PUTU,CORSRE, SUBTRC, AND  INTADD. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  /QDAT/NQSIZE , NQERR 
EXTERNAL  U1 ,F 
DIMENSION  EPS ( 10 ) , IR2 ( 10 ) 

COMMON  /SOLTAU/M , NGRSOL , PT 
C 
C 

C  SET  UP  ARRAYS  1  TO  M  FOR  THE  SOLUTIONS 

C  AND  ARRAYS  M+1  TO  2*M  FOR  THE  RIGHT  HAND  SIDES, 

C  AND  ARRAYS  2M+1  TO  3M-1  FOR  TAU  ARRAYS  AND 

C  SET  ASIDE  SPACE  FOR  GRID-7  SOLUTION  IN  3M=NGRSOL  GRID 

C  AND  CHECK  THAT  Q  ARRAY  IS  LARGE  ENOUGH 

NQERR=0 
DO  1  K=1 , M 
K2=2**(K-1 ) 

CALL  GRDFN (K,NX0*K2+1 ,NY0  *K2  +  1 , H0/K2 ) 

CALL  GRDFN (K+M,NX0*K2+1 ,NY0*K2+1 ,H0/K2 ) 

1  CALL  GRDFN (K+2  *M,NX0  *K2+1 ,NY0*K2+1 , H0/K2 ) 

NGRSOL=3*M 
PRINT  90,NQSIZE 

90  FORMAT ( '  SIZE  OF  Q  ARRAY  =  ',  110) 

IF ( NQERR. EQ.O)GOTO  92 
PRINT  91, NQERR 

91  FORMAT {'  ***  ERROR  IN  GRDFN  ***  ARRAY  Q  NOT  LARGE  ENOUGH  ***' 

*  /,'  ARRAY  Q  SIZE  SHOULD  BE  AT  LEAST  =' ,  110) 

STOP 

92  CONTINUE 
C 
C 

CALL  SOLRED 
C 

C  INITIALIZE 

WU=0 

CALL  PUTF( LIN , U1 , 0 ) 

DO  10  L=LIN,M 
C 

C  BEGIN  NEW  FINEST  LEVEL 

C 

PRINT  6,L 

6  FORMAT (1H0,60(1H. ) , 13 , 2X , 60 ( 1H. ) /) 

CALL  PUTF(L+M,F,2) 

TOLL=TOL* ( RATIO**L ) 

EPS { L ) =TOLL 

WU=.25*WU 

NCYCL=NCYC 

IF ( L . EQ . M ) NCYCL=NCYCM 

ICYC=0 

WMAXL=WMAX 

IF(L.EQ.M) WMAXL=WMAXM 
PRECL=PREC 

I F ( L . EQ . M) PRECL=PRECM 
C 
C 

K=L 

IR2 (L ) =0 
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328. 

329. 
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332. 
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C 

C 

C 


c 

c 


c 

c 


c 


c 


c 

c 


c 

c 


c 

c 

c 


BEGIN  A  NEW  WORK  LEVEL 


5  IR1=0 

ERR=1.E30 

RELAX  ONCE  ON  GRID  K 

3  ERRP=ERR 

CALL  RELAX (K, K+M, ERR) 

WU=WU+4 .** (K-L) 

IR1=IR1+1 

IR2(K)=IR2<K)+1 

WRITE (6,40 )K, ERR, WU, IR1 , IR2 (K ) 

40  FORMAT ( '  LEVEL', 12, ‘  RESIDUAL  NORM= ' ,  DIO. 3,'  WORK**' 
*  IR1=  ',12,'  IR2(K)=',I2> 

DECIDE  WHICH  GRID  TO  USE  NEXT 
IF  ( WU.GE. WMAXL )GOTO  20 
IF ( ERR. LT . EPS ( K )  ) GOTO  2 
IF( IR2 (K) .NE.NR2 )GOTO  8 
IF (  K.LT.L)GOTO  2 

ICYC=ICYC+1 

IF( ICYC.EQ.NCYCL  .AND.  L. NE.LIN )GOTO  20 
IF ( ICYC . EQ . NCYCLN  .AND.  L.EQ.LIN )GOTO  20 
IR2(L)=0 
IR1=0 

8  IF(IR1.EQ.NR1 )GOTO  4 

IF(IR1.EQ.1. OR.ERR.LT.  ERRP*ETA)GO  TO  3 

GO  TO  COARSER  GRID 

4  IF(K.EQ.1 )GOTO  3 

CALL  RESSW(K,K+M,K+M-1  ) 

CALL  RESBW(K,K+M,K+2*M-1) 

EPS(K-1 )=DELTA*ERR 
K=K-1 

CALL  PUTU ( K+ 1 , K ) 

CALL  CORSRE ( K, K+M) 

ITAUEX=0 

IF ( (ITAU.EQ. 1 ) .AND. (L.GT.LIN ) .AND. (K.EQ.L-1 ) )ITAUEX=1 
CALL  TAUC AM ( K , K+M , K+2  *  M , I TAUEX , TAUGNM) 

PRINT  60, TAUGNM, K 

60  FORMAT (50X, 'GREEN  NORM  OF  TAU-Z  = ' ,E12 . 3 , 5X , ’K=' , 12 ) 
IF(K.EQ. (L-1 ) ) EPS ( L ) =DMAX 1 ( PRECL  *TAUGNM , TOLL ) 

IR2 (K)=0 
GOTO  5 

GO  TO  FINER  GRID 
2  IF( K. EQ.L )GOTO  20 
CALL  SUBTRC( K+1 , K) 

CALL  INTSW(K,K+1 ) 

K=K+1 
GOTO  5 


FINISHED  WITH  LEVEL  L 
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343. 

344. 
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347. 
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349. 
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399. 


20  CONTINUE 
C 

C  THE  NEXT  SEVEN  STATEMENTS  COMPUTE  THE  GREEN  NORM  OF  TAU 

C  AND  THE  GREEN  AND  L-INFINITY  NORMS  OF  THE  ERROR 

C  (  IF  ACCURATE  SOLUTION  IS  KNOWN  ) 

11  CALL  RESSW(L,L+M,L+M-1 ) 

CALL  RESBW ( L , L+M , L+2  *M- 1  ) 

CALL  PUTU ( L,L-1 ) 

CALL  CORSRE(L-1 ,L-1+M) 

CALL  TAUCAM(L-1,L-1+M,L-1+2*M,0,TAUGNM) 

K=L-1 

PRINT  60,TAUGNM,K 
CALL  DIFFMX(L) 

C 

C 

IF(L.EQ.M)GOTO  10 
CALL  INTRP3 (L,L+1 ) 

CALL  PUTB( U1 ,L+1 ) 

C 

10  CONTINUE 
RETURN 
END 
C 
C 

SUBROUTINE  CORSRE (K,KRHS ) 

C  APPLIES  THE  DIFFERENCE  OPERATOR  ON  GRID  K 

C  TO  THE  GRID  FUNCTION  IN  ARRAY  K,  AND  ADDS  THE  RESULT  TO  THE 

C  VALUES  IN  ARRAY  KRHS. 

C  KRHS  KRHS  K  K,0 

C  B  =  R  +  A  U 

C 

C  THE  RESULT  IS  STORED  IN  ARRAY  KRHS. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  Q<  18000 ), 1ST < 200 ),IRHS( 200) 

CALL  KEY(K,IST,II,JJ,H) 

CALL  KEY (KRHS, IRHS, II, JJ,H) 

11=11-1 
J1=JJ-1 
DO  1  1=2,11 
IR=IRHS ( I ) 

IO=I ST(I) 

IM=IST( 1-1  ) 

IP=IST (1+1 ) 

DO  1  J=2,J1 

A=-Q { IR+J  )-Q ( IO+J+1 )-Q (IO+J-1 )-Q ( IM+J )-Q ( IP-KT) 

1  Q { IR+J )  =  -A-4  «*Q ( IO+J ) 

RETURN 

END 

C 

C 

SUBROUTINE  DIFFMX(K) 

C  NOT  TIMED 

C  COMPARES  SOLUTION  ON  GRID  K  WITH  ACCURATE  SOLUTION 

C  STORED  IN  GRID  NGRSOL 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  Q ( 18000 ) , 1ST (200 ) , ISTA(200 ) 

COMMON  /SOLTAU/M, NGRSOL , PT 
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=====  APPX-C-PFMG  ===== 

TIME=URTIMG ( 0 ) 

CALL  KEY (  K, 1ST, II , JJ ,  H  ) 

CALL  KEY (NGRSOL, ISTA, IIA, JJA, HA ) 

DIFMX=0 . 

DIFGNM=0 
SOLMX=0 . 

SOLGNM=0 

INTERV= ( IIA-1 )/( II-1  ) 

DO  1  1=1,11 
X=(I-1 )*H 
IA=(I-1)*INTERV+1 
DO  1  J=1,JJ 
Y=(J-1 ) *H 
JA=(J-1 )*INTERV+1 

DIF=ABS (Q ( ISTA( 1A )+JA)  -Q  ( 1ST ( I  )+J )  ) 

DIFGNM=DIFGNM+DIF*DIF 
SOL=ABS (Q ( ISTA( IA )+JA ) ) 

SOLGNM=SOLGNM+SOL*SOL 
SOLMX=AMAX 1 (SOL,SOLMX ) 

1  DIFMX=AMAX1 (DIF,DIFMX) 

DIFGNM=SQRT ( DXFGNM) /H 
PRINT  101 ,DIFMX,DIFGNM 

101  FORMAT ( 1 5X , '  SOLUTION  ERROR:  L  INFINITY  NORM  =',E13.5, 

*  5X , 1 GNORM  =  ' ,E13.5) 

SOLGNM=SQRT ( SOLGNM ) /H 
PRINT  1 02 ,SOLMX, SOLGNM 

102  FORMAT ( 1 5X , '  SOLUTION  :  L  INFINITY  NORM  =',E13.5, 

*  5X, 'GNORM  =  ' ,E13 . 5 ) 

PRINT  103, DIFMX/SOLMX , DIFGNM/ SOLGNM 

103  FORMAT ( 1 5X , '  RELATIVE  ERROR:  L  INFINITY  NORM  =',E13.5, 

*  5X, 'GNORM  =  ' ,E1 3.5) 

CALL  URT IMS (TIME ) 

RETURN 

END 

C 

C 

SUBROUTINE  GRDFN(N, IMAX , JMAX , HH ) 

C  SETS  UP  ARRAY  N. 

C  IMAX  THE  DIMENSION  IN  THE  X  DIRECTION 

C  JMAX  THE  DIMENSION  IN  THE  Y  DIRECTION 

C  HH  THE  GRID  SIZE 

C  THE  ARRAY  NST  CONTAINS  THE  STARTING  ADDRESSES  OF  THE  ARRAYS. 

C  THE  ARRAY  IMX  CONTAINS  THE  MAXIMUM  ROW  NUMBERS 

C  THE  ARRAY  JMX  CONTAINS  THE  MAXIMUM  COL  NUMBERS 

C  THE  ARRAY  H  CONTAINS  THE  GRID  SIZES. 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON /GRD/NST (20 ), IMX ( 20 ) , JMX ( 20 ) , H ( 20 } 

COMMON  /QDAT/NQSIZE , NQERR 
DATA  IQ/1/ 

NST (N )=IQ 
IMX (N ) =IMAX 
JMX (N)= JMAX 
H(N )=HH 

IQ=IQ+I MAX* JMAX 

IF( IQ. LE.NQSIZE+1 )  RETURN 

NQERR=IQ-1 
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END 

C 

c 

SUBROUTINE  INTSW(KC,KF) 

C  INTERPOLATES  CORRECTION  ON  COARSE  GRID  KC 

C  AND  ADDS  TO  SOLUTION  ON  GRID  KF. 

C  KF  KF  KC  KF  KF 

C  U  PHI (  I  W  +  U  i  U  ) 

C  KC 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  / SWDAT/NFGSW / NINTSW , NRESS W 
GOTO (1,2), NINTSW 
C 

1  CALL  INTADD ( KC ,  KF ) 

RETURN 

C 

2  CALL  INTADM ( KC ,  KF ) 

RETURN 

END 

C 

C 

SUBROUTINE  INTADD (KC , KF ) 

C  LINEARLY  INTERPOLATES  CORRECTION  ON  COARSE  GRID  KC 

C  AND  ADDS  TO  SOLUTION  ON  GRID  KF. 

C  KF  KF  KC  KF  KF 

C  U  PHI (  I  W  +  U  ;  U  ) 

C  KC 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  Q( 18000 ) ,ISTC(200),ISTF(200) 

CALL  KEY ( KC , ISTC , IIC , JJC ,  HC ) 

CALL  KEY ( KF, ISTF, IIF, JJF,HF) 

DO  1  IC=2 , IIC 

IF=2*IC-1 

JF=1 

IFO=ISTF(IF) 

IFM=ISTF( IF-1  ) 

ICO=ISTC(IC) 

ICM=ISTC(IC-1 ) 

DO  1  JC=2,JJC 
JF=JF+2 

A=. 5* (Q( ICO+JC) +Q ( ICO+JC-1 ) ) 

AM=. 5* (Q ( ICM+ JC ) +Q ( ICM+ JC- 1 ) ) 

Q(IFO+JF)  =  Q(IFO+JF)+Q( ICO+JC) 

Q ( IFM+JF)  =  Q ( IFM+JF )+. 5* (Q( ICO+JC )+Q ( ICM+JC ) ) 

Q( IFO+JF-1 )=Q ( IFO+JF-1 )+A 
1  Q( IFM+JF- 1 )  =  Q ( IFM+JF-1 )+. 5* ( A+AM) 

RETURN 

END 

C 

SUBROUTINE  INTADM ( KC ,  KF ) 

C  MODIFICATION  #6 . 

C  LINEARLY  INTERPOLATES  CORRECTION  ON  COARSE  GRID  KC 

C  AND  ADDS  TO  SOLUTION  ON  GRID  KF. 

C  CORRECTION  ONLY  ADDED  IF  SOLUTION  U  ON  FINE  GRID  IS 

C  NOT  ZERO.  SEE  (5.15). 
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C  KF  KF  KC  KF 

C  U  =  I  U  +  U 

C  KC 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  Q( 18000) , ISTC( 200 ) , ISTF(200 ) 

CALL  KEY  (KC,  ISTC ,  IIC ,  JJC ,  HC ) 

CALL  KEY ( KF, ISTF,IIF, JJF, HF) 

DO  1  IC=2 , IIC 

IF=2*IC-1 

JF=1 

IFO=ISTF( IF) 

IFM=ISTF ( IF-1  ) 

ICO=ISTC(IC) 

ICM=ISTC( IC-1 ) 

DO  1  JC=2 , JJC 
JF=JF+2 

A=. 5*(Q ( ICO+JC )+Q { ICO+JC-1 )) 

AM=. 5* (Q ( ICM+JC ) +Q ( ICM+JC-1 ) ) 

IF (Q( IFO+JF) .NE. 0 )Q{ IFO+JF)  =  Q (IFO+JF )+Q( ICO+JC) 

IF(Q ( IFM+JF) . NE . 0 )Q ( IFM+JF)  =  Q ( IFM+JF )+• 5* (Q ( ICO+JC )+Q ( ICM+JC ) ) 

IF ( Q ( IFO+JF-1 ) .NE. 0 )Q ( IFO+JF-1 )=Q(IFO+JF-1 )+A 

IF(Q ( IFM+JF- 1 ) . NE. 0 )Q  ( IFM+JF- 1 )  =  Q( IFM+JF- 1 )+.5*(A+AM) 

1  CONTINUE 
RETURN 
END 
C 
C 
C 

SUBROUTINE  INTRP3 (KC, KF) 

C  PERFORMS  CUBIC  INTERPOLATION 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  Q(18000) , IUF( 200 ) , IUC(200 ) 

CALL  KEY (KF, IUF , I IF, JJF, HF) 

CALI  KEY (KC,IUC,IIC,JJC,HC) 

C 

C  KF  KF  KC  KF 

C  U  =  J  U  +  U 

C  KC 

c 

C  INTERPOLATE  IN  COARSE  COLUMNS  USING  COARSE  COLUMN  DATA 

DO  20  IC=1,IIC 
IF=2*IC-1 
IFO=IUF( IF ) 

ICO=IUC(IC) 

Q ( IFO+1 ) =Q ( ICO+1 ) 

C  FIRST  POINT  IN  COLUMN.  USE  EQU  (6.3) 

Q (IFO+2)-(5*Q( ICO+1 )+15*Q(ICO+2)-5*Q(ICO+3)+Q(ICO+4) )/16 

JJC2=JJC-2 

DO  10  JC=2 , JJC2 

JF=2*JC-1 

Q ( IFO+JF ) «Q ( ICO+JC ) 

C 

C  INTERIOR  POINT  IN  COLUMN.  USE  EQU  (6.2) 

Q ( IFO+JF+ 1 )=(-Q( ICO+JC- 1 ) +9  *Q ( ICO+JC ) 

*  +9*Q(ICO+JC+1 )-Q(ICO+JC+2 ) )/16. 

1 0  CONTINUE 
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571. 

Q ( IFO+J JF-2 ) =Q ( ICO+JJC- 1 ) 

572. 

Q ( IFO+JJF-1 )=(Q( ICO+JJC-3 )-5  *Q ( ICO+JJC-2 ) 

573. 

C 

574. 

C 

LAST  POINT  IN  COLUMN.  USE  EQU  (6.3) 

575. 

★ 

+ 15*Q( ICO+JJC- 1 ) +5 *Q( ICO+JJC) )/16 

576. 

Q ( IFO+JJF ) =Q ( ICO+JJC ) 

577. 

20 

CONTINUE 

578. 

c 

579. 

c 

INTERPOLATE  IN  INTERMEDIATE  FINE  COLUMNS 

580. 

c 

USING  ROW  DATA 

581. 

c 

582. 

c 

FIRST  COLUMN.  USE  EQU  (6.3) 

583. 

IM1 =IUF( 1 ) 

584. 

IO=IUF ( 2 ) 

585. 

IP1=IUF( 3 ) 

586. 

IP3=IUF(5) 

587. 

IP5=IUF(7) 

588. 

DO  30  J=1,JJF 

589. 

Q(IO+J)=(5*Q(IM1+J)+15*Q(IPl+J) 

590. 

* 

-5*Q(IP3+J)+Q(IP5+J) )/16. 

591. 

30 

CONTINUE 

592. 

c 

593. 

c 

INTERMEDIATE  COLUMNS .  USE  EQU  (6.2) 

594. 

IIF3=IIF-3 

595. 

DO  40  I=4,IIF3,2 

596. 

IM3=IUF(I-3) 

597. 

IMl=IUF( 1—1 ) 

598. 

IO=IUF( I ) 

599. 

IPl=IUF( 1+1 ) 

600. 

IP3=IUF(I+3 ) 

601. 

DO  40  J=1 , JJF 

602. 

Q ( IO+J ) = ( -Q ( IM3+J ) +9*Q ( IM 1 +J ) 

603. 

* 

+9*Q(IP1+J)-Q(IP3+J))/16. 

604. 

40 

CONTINUE 

605. 

c 

606. 

c 

LAST  COLUMN.  USE  EQU  (6.3) 

607. 

IM5=IUF ( I IF-6 ) 

608. 

IM3=IUF( IIF-4 ) 

609. 

IMl=IUF(IIF-2) 

610. 

IO=IUF( IIF-1 ) 

611. 

IP1=IUF( IIF) 

612. 

DO  50  J=1  ,  JJF 

613. 

Q ( IO+J )=(Q( IM5+J ) -5  *Q ( IM3+J ) 

614. 

* 

+15*Q(IM1+J)+5*Q(IP1+J) )/16 

615. 

50 

CONTINUE 

616. 

RETURN 

617. 

END 

618. 

c 

619. 

c 

620. 

SUBROUTINE  KEY(K, 1ST, IMAX, JMAX,HH) 

621. 

c 

RECOVERS  THE  INFORMATION  ABOUT  ARRAY  K  SET  I 

622. 

c 

THE  SUBROUTINE  GRDFN. 

623. 

c 

THE  VALUE  OF  THE  GRID  FUNCTION  AT  THE  POINT 

624. 

c 

IS  ADDRESSED  AS  U(IST(J)+I). 

625. 

c 

626. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

627. 

COMMON /GRD/NST (20 ) , IMX (20 ) , JMX (20 ) , H ( 20 ) 

APPX-C-PFMG 


628. 

629. 

630. 

631. 

632. 

633. 

634. 

635. 

636. 
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638. 

639. 
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642. 

643. 

644. 

645. 

646. 

647. 

648. 

649. 

650. 

651. 

652. 

653. 

654. 

655. 

656. 

657. 

658. 

659. 

660. 
661. 
662. 

663. 

664. 

665. 

666. 

667. 

668. 

669. 

670. 
671  . 

672. 

673. 

674. 

675. 

676. 

677. 

678. 

679. 

680. 
681. 
682. 

683. 

684. 


DIMENSION  1ST ( 1  ) 

IMAX=IMX(K) 

JMAX=JMX { K ) 

IS=NST(K)-JMAX-1 
DO  1  1=1 ,IMAX 
IS=IS  +  JMAX 
1  1ST ( I )=IS 
HH=H ( K ) 

RETURN 

END 

C 

C 

SUBROUTINE  PUTB(F,K) 

C  INSERTS  THE  BOUNDARY  VALUES  OF  THE  FUNCTION  F 

C  EVALUATED  AT  THE  POINTS  OF  GRID  K 

C  INTO  THE  ARRAY  K. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  Q( 18000) 

DIMENSION  1ST (200) 

CALL  KEY ( K, 1ST, II , JJ , H ) 

111=11-1 
DO  1  J=1,JJ 
X=0. 

Y= ( J-1 ) *H 

Q(IST( 1 )+J )=F(X, Y) 

X=(II-1 )*H 

Q( 1ST (II )+J )=F(X, Y) 

1  CONTINUE 

DO  2  1=2,111 
Y=0. 

X=(I-1 )*H 
Q ( 1ST (I )+1 )=F(X,Y ) 

Y=( JJ-1 )*H 
Q ( 1ST (I )+JJ )=F( X, Y ) 

2  CONTINUE 
RETURN 
END 

C 

C 

C 

C 

SUBROUTINE  PUTF( K,F,NH ) 

C  INSERTS  THE  VALUES  OF  THE  FUNCTION  F 

C  EVALUATED  AT  THE  POINTS  OF  GRID  K 

C  AND  MULTIPLIED  BY  GRIDSIZE**NH 

C  INTO  THE  ARRAY  K. 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  Q ( 18000 ), 1ST (600) 

CALL  KEY  (K, 1ST, II , JJ , H) 

H2=H**NH 
DO  1  1=1,11 
DO  1  J=1,JJ 
X=(I-1 )*H 
Y=(J-1)*H 

1  Q(IST(I)+J)=F(X,Y)*H2 

RETURN 
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685. 

END 

686. 

C 

687. 

C 

688. 

SUBROUTINE  PUTU(KF,KC) 

689. 

c 

THIS  SUBROUTINE  INJECTS  THE  SOLUTION  ON 

THE 

FINE  GRID 

690. 

c 

KF  INTO  THE  COARSE  GRID  KC. 

691. 

c 

KC/0  KC  KF 

692. 

c 

U  =  I  U 

693. 

c 

KF 

694. 

c 

695. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

696. 

COMMON  Q( 18000 ) ,IUF(200 ) , IUC(200 ) 

697. 

CALL  KEY (KF, IUF , IIF, JJF , HF ) 

698. 

CALL  KEY (KC,IUC, IIC , JJC, HC ) 

699. 

DO  1  IC=1,IIC 

700. 

IF=2  *IC- 1 

701 . 

IFO=IUF( IF) 

702. 

ICO=IUC(IC) 

703. 

JF=- 1 

704. 

DO  1  JC=1,JJC 

705. 

JF=JF+2 

706. 

Q(  ICO+JC )  =  Q ( IFO+JF ) 

707. 

1  CONTINUE 

708. 

RETURN 

709. 

END 

710. 

c 

711. 

c 

712. 

SUBROUTINE  RELAX (K,KRHS, ERR) 

713. 

c 

NORMAL  RELAXATION 

714. 

c 

CARRIES  OUT  ONE  GAUSS-SEIDEL  PROJECTED 

715. 

c 

SWEEP  ON  THE  GRID  K  WITH  RIGHT  HAND  SIDE 

IN 

ARRAY  KRHS 

716. 

c 

RETURNS  WITH  ERR=  G-NORM  OF  THE  DYNAMIC 

RESIDUALS 

717. 

c 

718. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

719. 

COMMON  Q(18000) ,IST(200) ,IRHS(200) 

720. 

CALL  KEY (K, 1ST, II , JJ,H ) 

721. 

CALL  KEY (KRHS, IRHS , II , JJ,H ) 

722. 

11=11-1 

723. 

J1=JJ-1 

724. 

ERR=0 . 

725. 

DO  1  1=2,11 

726. 

IR=IRHS (I ) 

727. 

IO=IST (I ) 

728. 

IM=IST(I-1 ) 

729. 

IP=IST (1+1 ) 

730. 

DO  1  J=2,J1 

731. 

A=Q (IR+J )-Q(IO+J+1 )-Q ( IO+J-1 )-Q ( IM+J )-Q ( IP+J ) 

732. 

QT=-.25*A 

733. 

QN=MAX(0.0,QT) 

734. 

ERR=ERR+ (QN-Q ( IO+J ) ) * *2 

735. 

1  Q(IO+J)=QN 

736. 

ERR=SQRT ( ERR ) /H 

737. 

RETURN 

738. 

END 

739. 

c 

740. 

SUBROUTINE  RESBW ( KF , KRF , KRC ) 

741. 

c 

SAME  AS  RESSW  EXCEPT  THAT  ONLY  THE  RHS  B 

IS 

TREATED 
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742. 

743. 

744. 

745. 

746. 

747. 

748. 

749. 

750. 
751  . 

752. 

753. 

754. 

755. 

756. 

757. 

758. 

759. 

760. 

761. 

762. 

763. 

764. 

765. 

766. 

767. 

768. 

769. 

770. 
771  . 

772. 

773. 

774. 

775. 

776. 

777. 

778. 

779. 

780. 
781  . 

782. 

783. 

784. 

785. 

786. 

787. 

788. 

789. 

790. 

791. 

792. 

793. 

794. 

795. 

796. 

797. 

798. 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 


CALCULATES  THE  RESIDUAL  ON  GRID  KF  WITH  RIGHT  HAND  SIDE 
IN  ARRAY  KRF  ,  AND  TRANSFERS  INTO  ARRAY  KRC. 

BEFORE  TRANSFER,  THE  RESIDUAL  IS  SCALED 
BY  MULTIPLYING  BY  THE  FACTOR  4  TO  TAKE  ACCOUNT  OF  THE 
FACT  THAT  THE  GRID  SIZE  ON  GRID  KF  IS  HALF  THE 
GRIDS IZE  ON  GRID  KC. 

KRC  KC  KRF 

R  =  4*S  (  B  ) 

KF 

COMMON  /SWDAT/NFGSW,NINTSW,NRESSW 
GOTO  (1,2) ,NRESSW 

1  CALL  RESBAL ( KF , KRF , KRC ) 

RETURN 

2  CALL  RESBL 1 ( KF , KRF , KRC ) 

RETURN 

END 


SUBROUTINE  RESBAL ( KF , KRF , KRC ) 

SAME  AS  RESCAL  EXCEPT  THAT  ONLY  RHS  B  IS  TREATED 
CALCULATES  THE  RESIDUAL  ON  GRID  KF  WITH  RIGHT  HAND  SIDE 
IN  ARRAY  KRF  ,  AND  INJECTS  INTO  ARRAY  KRC. 

BEFORE  INJECTION,  THE  RESIDUAL  IS  SCALED 
BY  MULTIPLYING  BY  THE  FACTOR  4  TO  TAKE  ACCOUNT  OF  THE 
FACT  THAT  THE  GRID  SIZE  ON  GRID  KF  IS  HALF  THE 
GRIDSIZE  ON  GRID  KC. 

KRC  KC  KRF 

R  =  4  *S  (  B  ) 

KF 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  Q ( 18000 ) ,IUF( 200) ,IRF(200),IRC( 200) 

CALL  KEY ( KF , IUF , IIF , JJF , HF ) 

CALL  KEY ( KRF, IRF, IIF, JJF, HF) 

CALL  KEY (KRC , IRC, IIC, JJC , HC ) 

IIC1=IIC-1 
JJC1=JJC-1 
DO  1  IC=2 , IIC1 
ICR=IRC( IC ) 

IF=2*IC-1 

JF=1 

IFR=IRF( IF ) 

IFO=IUF( IF ) 

IFM=IUF( IF-1 ) 

IFP=IUF ( IF+1 ) 

DO  1  JC=2 , JJC1 
JF=JF+2 

1  Q ( ICR+ JC ) =4 . * ( Q ( IFR+JF ) ) 

RETURN 

END 


SUBROUTINE  RESBL 1 (KF, KRF, KRC ) 

SAME  AS  RESCL1  EXCEPT  THAT  ONLY  RHS  B  IS  TREATED 


J . 
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799. 

C 

MODIFICATION  #5  UPDATED  JUNE  23  1980 

800. 

C 

USES  WEIGHTED  RESIDUALS  NEAR  THE  BOUNDARY 

801. 

c 

CALCULATES  THE  RESIDUAL  ON  GRID  KF  WITH  RIGHT  HAND 

802. 

c 

IN  ARRAY  KRF  ,  AND  INJECTS  INTO  ARRAY  KRC. 

803. 

c 

BEFORE  INJECTION,  THE  RESIDUAL  IS  SCALED 

804. 

c 

BY  MULTIPLYING  BY  THE  FACTOR  4  TO  TAKE  ACCOUNT  OF 

805. 

c 

FACT  THAT  THE  GRID  SIZE  ON  GRID  KF  IS  HALF 

THE 

806. 

c 

GRIDSIZE  ON  GRID  KC. 

807. 

c 

KRC  KC  KRF 

808. 

c 

R  =  4*S  (  B  ) 

809. 

c 

KF 

810. 

c 

811  . 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

812. 

COMMON  Q ( 18000 ) , IUF(200 ) , IRF (200 ) ,IRC(200) 

813. 

DIMENSION  R(9) 

814. 

CALL  KEY ( KF, IUF, IIF, JJF , HF ) 

815. 

CALL  KEY (KRF, IRF, IIF, JJF, HF) 

816. 

CALL  KEY (KRC, IRC, IIC, JJC, HC ) 

817. 

IIC1=IIC-1 

818. 

JJC1=JJC-1 

819. 

DO  1  IC=2  ,  IIC  1 

820. 

ICR=IRC( IC ) 

821. 

IF=2  *IC-1 

822. 

JF=1 

823. 

IFR=IRF( IF) 

824. 

IFO=IUF( IF) 

825. 

IFM=IUF( IF-1  ) 

826. 

IFP=IUF( IF+1 ) 

827. 

DO  1  JC=2,JJC1 

828. 

JF=JF+2 

829. 

IF (Q ( IFO+JF) . EQ. 0 )GOTO  2 

830. 

IF(Q( IFP+JF+1 ) .GT.O  .AND.  Q ( IFP+JF-1 ) .GT. 0 

.AND. 

831. 

* 

Q(IFO+JF+1 ).GT.O  .AND.  Q( IFO+JF-1 ) .GT. 0 

.AND. 

832. 

* 

Q(IFM+JF+1 ) .GT.O  .AND.  Q ( IFM+JF-1 ) .GT. 0 

.AND. 

833. 

* 

Q ( IFM+JF  ) . GT • 0  .AND.  Q(IFP+JF  ).GT.O 

)GOTO  2 

834. 

N=0 

835. 

DO  3  11=1,3 

836. 

I=IF+I 1 -2 

837. 

DO  3  J1=1 ,3 

838. 

J=JF+J 1-2 

839. 

N=N+1 

840. 

IR=IRF ( I ) 

841. 

IO=IUF( I ) 

842. 

IM=IUF( 1-1 ) 

843. 

IP=IUF ( 1+1 ) 

844. 

S=Q ( IR+J ) 

845. 

IF(Q ( IO+J ) .EQ. 0 )S=0 

846. 

R(N)=S 

847. 

3 

CONTINUE 

848. 

Q( ICR+JC)=R( 5 )+ . 5* ( R ( 2 )+R(4 )+R(6 )+R(8)+ 

849. 

★ 

. 5  * ( R ( 1 ) +R ( 3 ) +R ( 7 ) +R ( 9 )  )) 

850. 

GOTO  1 

851. 

2 

Q ( ICR+JC ) =4 . *Q ( IFR+JF ) 

852. 

1 

CONTINUE 

853. 

RETURN 

854. 

END 
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SUBROUTINE  RESSW ( KF , KRF , KRC ) 

CALCULATES  THE  RESIDUAL  ON  GRID  KF  WITH  RIGHT  HAND  SIDE 
IN  ARRAY  KRF  ,  AND  TRANSFERS  INTO  ARRAY  KRC. 

BEFORE  TRANSFER,  THE  RESIDUAL  IS  SCALED 
BY  MULTIPLYING  BY  THE  FACTOR  4  TO  TAKE  ACCOUNT  OF  THE 
FACT  THAT  THE  GRID  SIZE  ON  GRID  KF  IS  HALF  THE 
GRIDSIZE  ON  GRID  KC. 

KRC  KC  KRF  KF  KF 

R  =  4*S  (  B  -  A  U  ) 

KF 

COMMON  /SWDAT/NFGSW , NINTSW , NRESSW 
GOTO  (1,2), NRESSW 

1  CALL  RESCAL (KF, KRF , KRC ) 

RETURN 

2  CALL  RESCL1 (KF, KRF, KRC) 

RETURN 


SUBROUTINE  RESCAL ( KF , KRF , KRC ) 

CALCULATES  THE  RESIDUAL  ON  GRID  KF  WITH  RIGHT  HAND  SIDE 
IN  ARRAY  KRF  ,  AND  INJECTS  INTO  ARRAY  KRC. 

BEFORE  INJECTION,  THE  RESIDUAL  IS  SCALED 
BY  MULTIPLYING  BY  THE  FACTOR  4  TO  TAKE  ACCOUNT  OF  THE 
FACT  THAT  THE  GRID  SIZE  ON  GRID  KF  IS  HALF  THE 
GRIDSIZE  ON  GRID  KC. 

KRC  KC  KRF  KF  KF 

R  =  4*S  (  B  -  A  U  ) 

KF 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z> 

COMMON  Q ( 1 8000 ) , IUF( 200 ) , IRF( 200 ) , IRC ( 200 ) 

CALL  KEY ( KF , IUF , IIF, JJF, HF) 

CALL  KEY ( KRF , IRF , IIF, JJF,HF) 

CALL  KEY ( KRC , IRC , I IC , JJC , HC ) 

IIC1=IIC-1 
JJC1=JJC-1 
DO  1  IC=2, IIC1 
ICR=IRC( IC) 

IF=2*IC-1 

JF=1 

IFR=IRF( IF ) 

IFO=IUF( IF) 

IFM=IUF( IF- 1 ) 

IFP=IUF( IF+1 ) 

DO  1  JC=2, JJC1 
JF=JF+2 

S=Q (IFO+JF+1 )+Q ( IFO+JF-1 )+Q ( IFM+JF)+Q ( IFP+JF ) 

Q ( ICR+JC ) =4 . * (Q ( IFR+JF ) -S+4 . *Q ( IFO+JF) ) 

RETURN 

END 
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913. 

914.  C 

915.  C 

916.  C 

917.  C 

918.  C 

919.  C 

920.  C 

921.  C 

922.  C 

923.  C 

924.  C 

925.  C 

926. 

927. 

928. 

929. 

930. 

931. 

932. 

933. 

934. 

935. 

936. 

937. 

938. 

939. 

940. 

941. 

942. 

943. 

944. 

945. 

946. 

947. 

948. 

949. 

950. 

951. 

952. 

953. 

954. 

955. 

956. 

957. 

958. 

959. 

960. 

961. 

962. 

963. 

964. 

965. 

966. 

967. 

968. 

969. 


SUBROUTINE  RESCLl ( KF , KRF , KRC ) 

MODIFICATION  #5  UPDATED  JUNE  23  1980 
USES  WEIGHTED  RESIDUALS  NEAR  THE  BOUNDARY 
CALCULATES  THE  RESIDUAL  ON  GRID  KF  WITH  RIGHT  HAND  SIDE 
IN  ARRAY  KRF  ,  AND  INJECTS  INTO  ARRAY  KRC. 

BEFORE  INJECTION,  THE  RESIDUAL  IS  SCALED 
BY  MULTIPLYING  BY  THE  FACTOR  4  TO  TAKE  ACCOUNT  OF  THE 
FACT  THAT  THE  GRID  SIZE  ON  GRID  KF  IS  HALF  THE 
GRIDSIZE  ON  GRID  KC. 

KRC  KC  KRF  KF  KF 

R  =  4*S  (  B  -  A  U  ) 

KF 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  Q ( 18000 ),IUF( 200 ),IRF( 200 ),IRC( 200) 

DIMENSION  R(9 ) 

CALL  KEY ( KF, IUF, IIF, JJF, HF ) 

CALL  KEY ( KRF , IRF , IIF , JJF , HF ) 

CALL  KEY ( KRC, IRC, IIC, JJC, HC ) 

IIC1=IIC-1 
JJC1=JJC-1 
DO  1  IC=2,IIC1 
ICR=IRC( IC ) 

IF=2*IC-1 

JF=1 

IFR=IRF( IF ) 

IFO=IUF( IF) 

IFM=IUF( IF-1 ) 

IFP=IUF( IF+1 ) 

DO  1  JC=2 , JJC1 
JF=JF+2 

IF(Q(IFO+JF) . EQ.O )GOTO  2 

IF(Q(IFP+JF+1 ) .GT.O  .AND.  Q(IFP+JF-1 ) .GT.O  .AND. 

*  Q( IFO+JF+1 ) .GT.O  .AND.  Q(IFO+JF-1 ) .GT.O  .AND. 

*  Q(IFM+JF+1 ) .GT.O  .AND.  Q ( IFM+JF-1 ) .GT.O  .AND. 

*  Q(IFM+JF  ) • GT . 0  .AND.  Q(IFP+JF  ).GT.O  )GOTO  2 

N=0 

DO  3  11=1,3 
I=IF+I1-2 
DO  3  J1=1 , 3 
J=JF+J1-2 
N=N+1 
IR=IRF(I ) 

IO=IUF( I ) 

IM=IUF(I-1 ) 

IP=IUF( 1+1 ) 

S=Q (IO+J+1 )+Q ( IO+J-1 )+Q ( IM+J )+Q ( IP+J ) 

S=Q ( IR+J ) +4  *Q { IO+J ) -S 
IF ( Q ( IO+J ) . EQ . 0 ) S=0 
R(N)=S 
3  CONTINUE 

Q(ICR+JC)=R(5)+.5*(R(2 )+R(4)+R(6)+R(8)+ 

1  . 5* (R{ 1 )+R( 3 ) +R(  7 )+R(  9 )  )) 

GOTO  1 

2  S=Q  ( IFO+JF+ 1  )+Q  ( IFO+JF-1  )+Q(IFM+JF)+Q(IFP+JF) 

Q ( ICR+JC ) =4 . * ( Q ( I FR+JF ) -S+4 . *Q ( IFO+ JF ) ) 

1  CONTINUE 
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970. 

RETURN 

971. 

END 

972. 

C 

973. 

C 

974. 

SUBROUTINE  SOLPRT(K,MPRINT) 

975. 

C 

NOT  TIMED 

976. 

C 

PRINTS  THE  ARRAY  K  ON  THE  SUBARRAY  MPRINT. 

977. 

C 

IF  K< MPRINT,  PRINTS  ENTIRE  ARRAY  K 

978. 

C 

979. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

980. 

COMMON  Q( 18000), QTEM( 100 ),IST(600) 

981  . 

TIME=-URTIMG(0  ) 

982. 

CALL  KEY  (MPRINT, 1ST, IIM,JJ,H) 

983. 

CALL  KEY  (K, 1ST, II, JJ, H ) 

984. 

INTERV=1 

985. 

IF ( K. GT. MPRINT ) INTERV=( II -1 )/<IIM-1 ) 

986. 

DO  20  J=JJ, 1,-INTERV 

987. 

L=0 

988. 

DO  10  1=1,11, INTERV 

989. 

c 

X  AND  Y  ARE  NOT  PRINTED  HERE,  BUT  ARE  COMPUTED 

990. 

c 

CASE  A  LATER  VERSION  NEEDS  THEM. 

991. 

X=(I-1 )*H 

992. 

Y= ( J-1 ) *H 

993. 

L=L+1 

994. 

QTEM(L ) =Q ( 1ST (I )+J ) 

995. 

10 

CONTINUE 

996. 

PRINT  * , (QTEM(LL) ,LL=1 ,L ) 

997. 

20 

CONTINUE 

998. 

CALL  URTIMS (TIME ) 

999. 

RETURN 

1000. 

END 

1001  . 

c 

1002. 

c 

1003. 

c 

1004. 

SUBROUTINE  SOLRED 

1005. 

c 

NOT  TIMED 

1006. 

c 

PUTS  ACCURATE  SOLUTION  INTO  GRID  NGRSOL 

1007. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

1008. 

COMMON  Q( 18000 ) , ISTA( 200 ) ,QTEM( 600 ) 

1009. 

COMMON  /SOLTAU/M, NGRSOL, PT 

1010. 

COMMON  /PRBDAT/Y1 ,Y2,A,R 

1011. 

COMMON  /SWDAT/NFGSW , NINTSW , NRESSW 

1012. 

TIME=URTIMG( 0 ) 

1013. 

CALL  KEY (NGRSOL, ISTA, IIA, JJA, HA) 

1014. 

c 

1015. 

GOTO ( 1,2) ,NFGSW 

1016. 

c 

1017. 

c 

DAM  PROBLEM 

1018. 

c 

ACCURATE  SOL  IS  DOUBLE  PRECISION  ON  GRID  M=7 

1019. 

c 

WITH  INITIAL  GRID  2X3 

1020. 

c 

STORED  IN  FILE  10. 

1021. 

1 

CONTINUE 

1022. 

MA=7 

1023. 

I IMA=2*  * ( MA- 1 ) *2+1 

1024. 

INTERV=(IIMA-1 ) /( IIA-1 ) 

1025. 

JJMA= ( JJA- 1 ) *INTERV+ 1 

1026. 

REWIND  10 

94 
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1027. 

DO  20  JA=1 ,  JJMA 

1028. 

READ (10)  (QTEM(IA) ,IA=1,IIMA) 

1029. 

J=( JA-1 ) /INTERV+1 

1030. 

IF(  (J-1 )*INTERV  .NE.  JA-1  )GOTO  20 

1031. 

DO  10  1=1, IXA 

1032. 

IA=(I-1 )*INTERV+1 

1033. 

Q (ISTA(I )+J )=QTEM( IA ) 

1034. 

10 

CONTINUE 

1035. 

20 

CONTINUE 

1036. 

GOTO  1000 

1037. 

C 

1038. 

C 

1039. 

C 

PROBLEM  OF  SECTION  5:  (5.3)  AND  (5.4) 

1040. 

C 

EXACT  SOLUTION  KNOWN 

1041. 

2 

CONTINUE 

1042. 

D=2.5*R 

1043. 

DO  30  1=1, IXA 

1044. 

IO=ISTA( I ) 

1045. 

DO  25  J=1 , JJA 

1046. 

X=(I-1 )*HA 

1047. 

Y=( J-1 )*HA 

1048. 

A=DMAX1 ( 0 .DO ,D-R*X-Y ) 

1049. 

B=X+Y 

1050. 

G=A* A* ( DCOS ( B ) +2 ) 

1051. 

Q ( IO+J )=G 

1052. 

25 

CONTINUE 

1053. 

30 

CONTINUE 

1054. 

GOTO  1000 

1055. 

C 

1056. 

1000 

CONTINUE 

1057. 

CALL  URTIMS (TIME ) 

1058. 

RETURN 

1059. 

END 

1060. 

C 

1061. 

SUBROUTINE  SUBTRC ( KF , KC ) 

1062. 

C 

THIS  SUBROUTINE  COMPUTES  THE  VALUE  INJECTED  FROM  GRID  KF  TO 

1063. 

C 

GRID  KC  AND  SUBTRACTS  IT  FROM  THE  SOLUTION  ON  GRID  KC. 

1064. 

C 

KC  KC  KC  KF 

1065. 

C 

W  =  U  -  I  U 

1066. 

c 

KF 

1067. 

c 

1068. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

1069. 

COMMON  Q ( 18000 ) ,IUF(200 ) , IUC( 200 ) 

1070. 

CALL  KEY ( KF , IUF , IIF , JJF , HF ) 

1071. 

CALL  KEY ( KC, IUC , IIC, JJC,HC) 

1072. 

DO  1  IC=1,IIC 

1073. 

IF=2*IC-1 

1074. 

IFO=IUF( IF ) 

1075. 

ICO=IUC(IC) 

1076. 

JF=- 1 

1077. 

DO  1  JC=1,JJC 

1078. 

JF=JF+2 

1079. 

Q ( ICO+JC )  =Q  ( ICO+ JC ) -Q ( IFO+JF ) 

1080. 

1 

CONTINUE 

1081. 

RETURN 

1082. 

END 

1083. 

c 
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1084. 

C 

1085. 

SUBROUTINE  TAUCAM ( KU  ,  KR ,  KF , ITAU , TAUGNM ) 

1086. 

C 

COMPUTES  TAU  AND  TAU-Z  GREEN  NORM 

1087. 

C 

UPDATED  AUGUST  26  1980 

1088. 

c 

PERFORMS  TAU  EXTRAPOLATION  IF  ITAU-1 

1089. 

c 

BY  ADDING  TAU  TO  RHS  ON  GRID 

1090. 

c 

GRID  KU  CONTAINS  U 

1091 . 

c 

GRID  KR  CONTAINS  SUM  OF  FIRST  TWO  TERMS  IN  (6.7) 

1092. 

c 

PREVIOUSLY  OBTAINED  USING  RESSW  AND  CORSRE 

1093. 

c 

GRID  KF  CONTAINS  THIRD  BRACKET  IN  (6.7)  PREVIOUSLY 

1094. 

c 

COMPUTED  BY  RESBW 

1095. 

c 

ITAU  IS  PARAMETER  WHICH  DETERMINES  WHETHER  EXTRAPOLATION 

1096. 

c 

WILL  BE  PERFORMED 

1097. 

c 

TAUGNM  IS  RETURNED  AS  GREEN  NORM  OF  TAU-Z 

1098. 

c 

1099. 

c 

K-1  PT  K-1  K  K  K  K-1  K-1  K  K-1  1 

1100. 

c 

T  =2  *  (4S  (B  -  A  U  )  )  +  (A  I  U)  -  (4S  B 

1101. 

c 

-  K  K  K 

1102. 

c 

2**PT- 1 

1103. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

1104. 

COMMON  Q(18000),IKR(200),IKF(200),IKU(200) 

1105. 

COMMON  /SOLTAU/M , NGRSOL , PT 

1106. 

CALL  KEY(KR,IKR,II,JJ,HK) 

1107. 

CALL  KEY(KF,IKF,II,JJ,HK) 

1108. 

CALL  KEY (KU, IKU, II » JJ,HK) 

1109. 

A=2.**PT/(2.**PT-1 ) 

1110. 

TAUGNM=0 

1111. 

111=11-1 

1112. 

1113. 

DO  1  IK=2 ,111 

1114. 

IRK=IKR( IK) 

1115. 

IFKO=IKF ( IK ) 

1116. 

IO=IKU ( IK) 

1117. 

IM=IKU (IK-1 ) 

1118. 

IP=IKU ( IK+1 ) 

1119. 

DO  1  JK=2  ,  JJ 1 

1120. 

T=Q ( IRK+JK ) -Q ( IFKO+JK ) 

1121. 

T=A*T 

1122. 

IF(Q ( IO+JK) .EQ.O )T=0 

1 123. 

TAUGNM=TAUGNM+T*T 

1124. 

IF(  Q(JK+IO+1 ) .EQ.O  .OR. 

1125. 

*  Q ( IO+JK- 1 ) . EQ . 0  .OR.  Q ( IM+JK) .EQ. 0  .OR. 

1126. 

*  Q ( IP+JK) .EQ.O )  T=0 

1127. 

IF ( ITAU . EQ . 1 )Q ( I RK+JK )=T+Q ( IFKO+JK) 

1128. 

1  CONTINUE 

1129. 

TAUGNM=SQRT { TAUGNM ) /HK 

1130. 

RETURN 

1131. 

END 
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